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Abstract 

The impact of interparticle correlations on the behavior of Bose-Einstein Condensates 

(BECs) is discussed using two approaches. In the first approach, the wavefunction of 

a BEC is encoded in the iV-particle sector of an extended “catalytic state”. Going 

to a time-dependent interaction picture, we can organize the effective Hamiltonian 
_ 1 _ 1 

by powers of N 2 . Requiring the terms of order N 2 to vanish, we get the Gross- 
Pitaevskii Equation. Going to the next order, N °, we obtain the number-conserving 
Bogoliubov approximation. Our approach allows one to stay in the Schrodinger 
picture and to apply many techniques from quantum optics. Moreover, it is easier to 
track different orders in the Hamiltonian and to generalize to the multi-component 
case. In the second approach, I consider a state of N = l x n bosons that is 
derived by symmetrizing the n-fold tensor product of an arbitrary /-boson state. 
Particularly, we are interested in the pure state case for 1 = 2, which we call the 
Pair-Correlated State (PCS). I show that PCS reproduces the number-conserving 
Bogoliubov approximation; moreover, it also works in the strong interaction regime 
where the Bogoliubov approximation fails. For the two-site Bose-Hubbard model, I 
find numerically that the error (measured by trace distance of the two-particle RDMs) 



of PCS is less than two percent over the entire parameter space, thus making PCS a 
bridge between the superfluid and Mott insulating phases. Amazingly, the error of 
PCS does not increase, in the time-dependent case, as the system evolves for longer 
times. I derive both time-dependent and -independent equations for the ground state 
and the time evolution of the PCS ansatz. The time complexity of simulating PCS 
does not depend on N and is linear in the number of orbitals in use. Compared to 
other methods, e.g, the Jastrow wavefunction, the Gutzwiller wavefunction, and the 
multi-configurational time-dependent Hartree method, our approach does not require 
quantum Monte Carlo nor demanding computational power. 
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Chapter 1 
Introduction 


Today, if you have a demanding job for light, you use an optical laser. In 
the future, if there is a demanding job for atoms, you may be able to use 
an atom laser. 

- Wolfgang Ketterle 

When Charles H. Townes, Nikolay Basov, Aleksandr Prokhorov, and Theodore 
Maiman were working on masers and lasers in the 1950s and 60s, they probably did 
not anticipate the extremely wide range of applications of their work. Lasers are 
ubiquitous nowadays: from supermarket scanners to remote sensing, from CD players 
to holographic technologies. But what makes lasers so useful besides their extremely 
high powers? Compared to previous light sources, lasers are both monochromatic and 
coherent; these properties enable one to take advantage of—much more efficiently— 
the interference effects of light. Interference is the key to manipulating waves just- 
like Newton’s laws of motion are crucial to controlling the movements of particles. 
Lasers allow one to control the waveforms of light. 

Perhaps one of the most- surprising discoveries of quantum mechanics is the 
wave-particle duality, or equivalently, that any quantum particle can interfere with 
itself. This “weird” idea of de Broglie turned out- to be quite useful beyond its 
purely theoretical interest; e.g., the diffraction patterns of scattered neutrons can 
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provide information about the positions and movements of the nuclei in a sample. 
“Atomic physics after about 1975 has been interested in controlling the atom’s external 
quantum state ... . The ultimate payoff of external state manipulation is atom 
interferometry which involves making and reading out superpositions of the external 
quantum state,” said David E. Pritchard [Pri02], Atom interferometers [Kas02, 
CSP09], since their first demonstration [CM91], have had an impact on many fields of 
science and engineering, for example, acceleration and rotation sensors, gravitational- 
wave detectors, and various other force and field sensors. Despite their relatively 
short history, atomic sensors can already compete with the best classical or laser- 
based sensors [GusOO, KKD11], A typical principle of operation is as follows: A 
thermal cesium atomic beam crosses three laser interaction regions where two-photon 
stimulated Raman transitions between cesium ground states transfer momentum to 
atoms and divide, deflect, and recombine the atomic wavepackets; rotation induces a 
phase shift between the two possible trajectories and causes a change in the detected 
number of atoms with a particular internal state. 

The advantages of atom interferometers include short de Broglie wavelength, 
narrow frequency response, and long interrogation times. The trapping of atoms also 
allows the possibility of steadily splitting the wavefunction in coordinate space, for 
example the double-well configuration in Fig. 1.3; this property, which is unprecedented 
with light interferometers, is essential to experimental studies of spatially varying 
fields such as gravitational forces. 


1.1 Bose-Einstein Condensates 

As in the case of light, one needs a unified army of atoms, all marching in step—i.e., 
an atom laser—to enhance the performance of a matter-wave interferometer, instead 
of just an ensemble of uncorrelated atoms. An atom laser can be thought as a 
monochromatic and coherent matter-wave beam [BHE99]; monochromatic means 
that most of the atoms in the beam occupy the same quantum state, and coherent 
means that the phases of atom beams are well defined and can be correlated. Lasers 
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are made by stimulated emission, which is not an option for atoms; so how can 
we make an atom laser? Remember that a laser, where photons share the same 
quantum state, is a state of light with extremely low entropy; similarly, the entropy 
of an atom laser must be very low. One way to get low-entropy states is to cool 
the atoms, and the temperature at which an atom beam becomes an atom laser is 
called the critical temperature T c . It turns out that for classical (distinguishable) 
particles, the critical temperature T c goes to zero when the number of particles 
N is very large (the thermodynamic limit). Microscopic particles are not classical 
particles; instead they are indistinguishable particles obeying either Bose-Einstein 
statistics (bosons) or Fermi-Dirac statistics (fermions). Bosons tend to occupy the 
same quantum state, while fermions can only occupy different quantum states. In 
1925, based on the earlier work of Satyendra Nath Bose, Albert Einstein pointed out 
that a large fraction of particles in a boson gas condense to the lowest quantum state 
at a finite critical temperature T c , which depends only on the mass of the particles 
and the density of particles at the thermodynamic limit. This phenomenon, called 
a Bose-Einstein Condensate (BEC), is a consequence of the Bose-Einstein statistics 
rather than attractive interactions between the particles; BECs happen when the de 
Broglie waves of atoms overlap, and this gives an estimate of the critical temperature 



( 1 . 1 ) 


where h is the reduced Planck constant, k B is the Boltzmann constant, m is the mass 
per particle, N is the number of particles, and V is the volume of the gas. 

To get the ultimate resolution of an atom interferometer, we need a single-mode, 
coherent source that is also very bright, i.e., a Bose-Einstein condensate. Yet BECs 
are so fragile that they have not been observed in nature. Erwin Schrodinger 
wrote [Sch52], “The densities are so high and the temperatures so low ... the van 
der Waals corrections are bound to coalesce with the possible effects of degeneration 
.. . .” What Schrodinger meant was that the atoms will condense into a solid or a 
liquid before they become a BEC, but what he forgot to consider was a dilute gas 
in a metastable phase; the challenge is to find a window of the critical temperature 
T c accessible by cryogenics, while making the lifetime of the metastable state long 
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enough. It turned out that the feasible critical temperature T c is around or below 
one microkelvin; e.g., the first BEC realized at the NIST-JILA lab by Eric Cornell 
and Carl Wieman [CW02] had a temperature of 170 nK. 
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Figure 1.1: Temperature gradient from a BEC to the center of the sun. 


Conventional cryogenics such as the dilution refrigerators provide cooling to the 
temperature of superfluid 3 He, which is about 2 mK, but this temperature is still 
about 10, 000 times hotter than a BEC. In addition, it takes a long time to reach these 
very low temperature using conventional cryogenics, not to mention the associated 
high maintenance costs. To chill the atoms more efficiently, physicists invented laser 
cooling, which soon became a central topic of modern cryogenics. But how can lasers, 
emitting high-energy photons, be used for cooling? The answer again lies in the low 
entropy of lasers; even though they emit photons whose energy is about the same as 
those emitted by the Sun, the laser photons are actually very cold and thus able to 
extract entropy from the atoms as the photons heat up. The most common method 
of laser cooling, Doppler cooling, works via the recoil forces exerted by the light on 
the atoms; no matter which way an atom drifts, it always runs into a laser beam that 
slows it down. But Doppler cooling has its limit due to the natural linewidth of the 
atoms in use, which for Rubidium 85 is around 150 pK. Other more sophisticated 
laser cooling techniques, such as Sisyphus cooling and polarization gradient cooling, 
are able to approach the much lower recoil temperature; the recoil temperature, which 
is usually on the order of 1 pK, equals the recoil energy deposited in a single atom 
initially at rest by the spontaneous emission of a single photon. The final step that 
takes laser cooled atoms to BECs is evaporative cooling, which works by allowing 
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the atoms with highest energy to escape from the trap where the atoms are stored. 
There are many excellent reviews of the theory and experiment of dilute-gas atomic 
BECs [Tau94, DGPS99, TKS97, Ket02, CW02], 

1.2 Interferometers with BECs 

Physics with coherent matter waves is an emerging research held [Mpl03, BS04], 
The first observation of interference between two BECs was made in Ketterle’s 
group [ATM + 97], owing to the high density of atoms in their Ioffe-Pritchard trap. 
Anderson et al. [AK98] observed interference of BECs tunneling from an array of traps. 
Stenger et al. [SIC + 99] developed Bragg spectroscopy to measure properties of a 
condensate for high-resolution velocimetry. Shin et al. [SSP + 04] succeeded in making a 
BEC based interferometer by continuously deforming between single- and double-well 
trapping potentials; the device was then miniaturized by using an atom chip [SSJ + 05, 
SHA + 05]. 1 Hadzibabic et al. [HSB + 04] observed matter-wave interference between 
30 BECs with uncorrelated phases. Recently, a matter-wave interferometer that 
uses optical ionization gratings has also been demonstrated [HDG + 13]; this does not 
require any particular internal level structure and is thus universally applicable to 
different kinds of particles. In addition, the first controllable atomic circuit that 
functions analogously to a Superconducting Quantum Interference Device (SQUID) 
has been implemented [WBL + 13]. 

Before diving into atom interferometers, we review the five essential steps of a 
general interferometer: (i) prepare the initial state; (ii) split the initial state into a 
coherent superposition of two states; (iii) apply interactions that affect the two states 
differently, for example, due to their different spatial locations; (iv) recombine the 
evolved states coherently; and (v) measure the shift of the interference fringes. In the 

1 See [HVB01] for a theoretical proposal to use atom chips as a substitute for magneto¬ 
optical traps: “A versatile miniature de Broglie waveguide is formed by two parallel 
current-carrying wires in the presence of a uniform bias field ... it offers a remarkable range 
of possibilities for atom manipulation ... include controlled and coherent splitting of the 
wave function as well as cooling, trapping, and guiding.” 
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following, I briefly review two different configurations of atom interferometers based 
on the above prescription: the atoms are not required to be condensed into a BEC 
for the first configuration, while BECs are required for the second configuration. 


In the first configuration, moving clouds of ultracold atoms are subjected to three 
laser beams as in Fig. 1.2. The lasers are tuned to induce Raman transitions between 

tt/2 tt tt/2 



t = 0 t — T t = 2T + 5 


Figure 1.2: A typical atom interferometer: The first 7t/2 pulse transfers half of 
the atoms coherently from the ground state | g ) into the other hyperhne state 
| e). After flying apart for time T, the atom states are mirrored by a 7r pulse. 
The atoms overlap again at time t = 2T, and a second tt/2 pulse recombines the 
two clouds. The number of atoms left in the state | g ) depends on the phase 
difference accumulated between the two paths. Furthermore, if the atoms are 
cold enough to form a BEC and the time delay 6 is nonzero, one should be able 
to observe interference fringes for both of the two clouds at the detectors. 


the two hyperhne ground states of atoms. The first n/2 pulse transfers half of the 
atoms coherently from the ground state j g) into the other hyperhne state | e); the 
atoms excited to \e) receive a recoil force and separate from those still left in the 
ground state | g). After evolving for a time T, the atom states are mirrored by a n 
pulse, and the recoil force swaps the velocity of the atoms in the two hyperhne states. 
The atoms overlap again at time t = 2T, and a second tt/2 pulse recombines the two 
clouds. The number of atoms left in the state | g ) depends on the phase difference 
accumulated between the two paths. Furthermore, if the atoms are cold enough to 
form a BEC and the time delay 5 is nonzero, one should observe interference fringes 
for both of the two clouds at the detectors. 
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In the other configuration [SSP + 04], where the atoms are put in a double-well 
trapping potential, a BEC is considered necessary for the fringe patterns to be visible. 
The five steps to perform this interferometer are listed in Fig. 1.3: (i) cool the atoms 
trapped in a single-well potential to form a BEC; (ii) split the condensate by slowly 
deforming the single-well potential to a double-well potential; (iii) apply ac Stark shift 
potentials to either of the two separated condensates; (iv) turn off the double-well 
trapping potential, thus letting the condensates ballistically expand, overlap, and 
interfere; and (v) take an absorption image. In this configuration, the atoms are 


(i) (“) 



Figure 1.3: A double-well atom interferometer [SSP + 04]: (i) cool the atoms 
trapped in a single-well potential to form a BEC; (ii) split the condensate by 
slowly deforming the single-well potential to a double-well potential; (iii) apply 
ac Stark shift potentials to either of the two separated condensates; (iv) turn off 
the double-well trapping potential, and let the condensates ballistically expand, 
overlap, and interfere; and (v) take an absorption image. 

confined in a trapping potential until the measurement is made, while the atoms 
are free to fly in the first configuration. This confinement of atoms allows long 
interrogation time, and the relative phase of two condensates can be measured using 
a small sample of the atoms. In addition, confined atom interferometers, especially 
those using atom chips, can be small and portable. On the other hand, confined 
atom interferometers usually operate with high density to achieve large signals, 
from which several disadvantages follow: (i) Three-body interactions can cause the 
atomic gas to form a liquid or solid; (ii) strong two-body interactions can cause large 
depletion of the condensate, even at zero temperature; (iii) the matter-wave dynamics 
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becomes nonlinear; (iv) interactions suppress the Josephson oscillation and cause 
phase diffusion; (v) the potential wells have to be controlled very accurately in stiffness 
and depth to prevent additional frequency shifts; (vi) the collective excitations of 
the condensate have to be controlled carefully, because, for example, sound or shape 
oscillations may arise if the potential changes too suddenly. 


1.3 Nonclassical Aspects of BECs 

To describe the interference of BECs, perhaps the most obvious choice is to modify 
the existing theory for laser light fields. There are, however, two intrinsic difference 
between lasers and BECs: (i) Lasers are superpositions of different number states 
of photons, while the number of particles in a BEC is usually fixed; (ii) lasers are 
composed of noninteracting photons, while particles in a BEC generally interact with 
one another. In this section, I discuss how these two differences are handled in the 
literature. 

Similar to a laser light, a BEC is usually described by a coherent state, i.e., an 
eigenstate of the annihilation operator for a particular state of the single-atom Hilbert 
space. A well defined amplitude a and phase angle 0 = arg(a) is associated with this 
coherent state. When particle loss is negligible, however, the real condensate is much 
closer to a number state than to a coherent state. A number state is distributed evenly 
across all phase-plane directions, and no definite phase can be attributed to it. The 
fictitious phase to the coherent state breaks the rotational symmetry and introduces 
Goldstone bosons into the Bogoliubov approximation to the condensate [LY96] ; these 
should be treated as a defect of the mathematical description and not as a physical 
property of the atomic system. The Goldstone mode causes the condensate state 
to deviate linearly in time from a single condensate in a coherent state (i.e., this 
is a secular deviation, not an oscillation). This problem is particularly pesky when 
the condensate is in a trapping potential, where the Goldstone mode is a mixture of 
the condensate mode and modes orthogonal to it and thus cannot be removed easily. 
The solution to getting rid of the unphysical Goldstone mode is to adhere to the 
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fact that the condensate has a fixed number of particles, i.e., by using a Bogoliubov 
approximation where particle number is conserved [GA59, Gar97, CD98]. 

Due to the interparticle interactions, a BEC based interferometer is very different 
from a laser based interferometer. A nonlinear effect known as macroscopic quan¬ 
tum self-trapping in bosonic Josephson junctions has been discussed and observed 
in [SFGS97, SROO, CBF + 01, AGF + 05, AB06]. The collapse and revival of interference 
patterns of matter waves were observed in optic lattices [GMHB02, SDZ11]. It is 
worth mentioning that sometimes the nonlinear effect is useful to achieving high sensi¬ 
tivity. One example is that of fast moving solitons, created in the recombination stage, 
which can enhance the sensitivity of phase measurements [NH04], Nonlinear effects 
in BECs can be useful for parameter estimation problems [CS08, BDD + 09, TBD + 10]. 
In particular, they can be used to produce so-called spin squeezed states, which are 
useful for overcoming the shot-noise limit [Cav81, BK97, OTF + Ol, DBB02, GLM04, 
EGW+08, PS09, GZN+10, RBL+10, Grol2], 

To quantify the impact of interparticle interactions on a BEC based interferometer, 
many authors have used the two-mode model. Milburn et al. [MCWW97] showed 
that the mean-held solution is modulated by a quantum collapse and revival sequence. 
In [CLMZ98, SC98, SDCZ01], the authors argued that it is feasible to prepare, 
control, and detect a Schrodinger cat state with a two-component. BEC. Spekkens 
and Sipe [SS99] considered the transition from a single condensate to a fragmented 
condensate as the central barrier in a double-well trapping potential is raised. Menotti 
et al. [MACZ01] explicitly considered the spatial dependence of the mode functions. 
Mahmud et al. [MPR05] analyzed the two-mode model in the quantum phase-space 
picture. Several authors [HM08, LFS12] discussed how to optimally create quantum 
superpositions and squeezing using the bosonic Josephson Hamiltonian. 2 More 
realistically, the influences of the noncondensed modes on the relative phase of the 
two condensed modes were considered in [VL99, SCW09, GGMBS14], 

2 

The ground and excited states of the bosonic Josephson Hamiltonian can be solved 
exactly by using the algebraic Bethe ansatz [LZMG03, ZLMG03]. 
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Interesting phenomena such as quantum phase transitions from a superfluid to a 
Mott insulator in a Bose gas [JBC + 98, GME + 02] cannot be interpreted within the 
mean-field approach. Also, a proper analysis of the atom interferometers discussed 
above requires one to go beyond mean-held theory. The most common way to do 
so is by perturbing the mean held with collective excitations, i.e., the Bogoliubov 
approximation. When interparticle interactions are strong or the mean held is unstable, 
the Bogoliubov approximation fails. In this section, I introduce some nonperturbative 
methods that go beyond the mean held theory. 

One powerful idea to deal with strong interactions is to construct a many-particle 
wavefunction from a two-particle state. In 1950s, Jastrow introduced the wavefunction 
which bears his name; the IV-particle state is a product of two-particle states of all 
N(N — l)/2 pairs, 

N 

^(xi,x 2 ,...,x iV ) ~ JJ/(x i -x jfc ). (1-2) 

j,k =1 
j<k 

Many famous wavefunctions are of Jastrow type, e.g., the Laughlin wavefunc¬ 
tion [Lau83] and the Gutzwiller wavefunction [Gut63, RK91]. The Jastrow wave- 
function has found wide application in strongly interacting systems: It is used in 
variational quantum Monte Carlo as a trial wavefunction [UWW88]; it is used to 
show that the one-particle reduced density matrix of 4 He is an extensive quantity, 
thus verifying that BEC underlies superfluidity [Rea69]; and it is used to investigate 
the effect of interatomic correlations and the accuracy of the Gross-Pitaevskii equa¬ 
tion [FP99, DG01, CHM + 02], The validity of the Jastrow wavefunction is discussed 
in [KKLZ91], where the authors constructed a class of interacting-boson Hamiltonians 
whose exact ground-state wavefunctions are of Jastrow form. 

Another ansatz that many authors have adapted to discuss fragmentation of BECs 
is the Double-Fock State (DFS), or sometimes Twin-Fock State (TFS); most generally, 
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one can use a many-Fock state that takes the form 



where v is the number of fragments. Using this ansatz, Streltsov et al. [SCM04] 
argued that fragmentation of the ground state of a BEC only happens when the total 
number of particles is finite; Mueller et al. [MHUB06] showed that as degeneracies 
multiply, so do the varieties of fragmentation; and Alon et al. [ASC07] generalized 
the Gross-Pitaevskii equation to include multiple orbitals. Other authors treated 
fragmented BECs by evolving the single-particle reduced density matrix with the 
Bogoliubov back-reaction approximation [VA01, TAV07]. 

One interesting numerical method to simulate the dynamics of N interacting 
bosons is the Multiconhgurational Time-Dependent Hart.ree Method for Bosons 
(MCTDHB) [MMR05, ASC08]; instead of fixing the single-particle orbitals, MCTDHB 
uses time-dependent orbitals, which reduces the number of orbitals required to reach 
a certain precision. However, it still has to simulate the full quantum dynamics in a 
D-dimensional Hilbert space, with D = — N u_1 , where v is the number of 

orbitals. 


1.5 Outline of This Dissertation 

Chapter 2 is an introduction to the basics of BECs. The Gross-Pitaevskii (GP) 
equation is derived using an approach that “projects” the evolved state back to the 
product manifold; I show that while the error introduced by the “projection” is large 
for the entire wavefunction, it is quite small for the few-particle Reduced Density 
Matrices (RDMs). This explains why the GP equation works so well in practice, 
because only the few-particle RDMs—not the whole wavefunction—can be measured 
in a lab. Also, the relative phase of two BECs is discussed, with a brief explanation 
of the MIT experiment [ATM + 97]. 

In Chap. 3, I introduce an alternative way to derive the number-conserving 
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Bogoliubov approximation, where the many-body wave function of a BEC is “encoded” 
in the iV-particle sector of an extended catalytic state , which is a coherent state for 
the condensate mode and an arbitrary state for the noncondensed modes. By going to 
a time-dependent interaction picture, the coherent state is displaced to the vacuum, 
where all the held operators are small compared to N 1 ^ 2 . The resulting Hamiltonian 
can then be organized by powers of N~ L ^ 2 . Requiring the terms of order N 1 ^ 2 to 
vanish, we get the Gross Pitaevskii equation for the condensate wave function. Going 
to the next order, N°, we are able to derive equations equivalent to those found 
by Castin and Dum [CD98] for a number-conserving Bogoliubov approximation. In 
contrast to other approaches, the extended-catalytic-state approach allows us to 
calculate the state evolution in the Schrodinger picture instead of the Heisenberg 
picture. In addition, many techniques from quantum optics, such as quasi-probability 
distributions, can be used directly in our approach. Moreover, it is much easier to 
track different orders in the Hamiltonian with our approach, which allows one to go to 
approximations beyond second order. Last but not the least, the number-conserving 
Bogoliubov approximations for multicomponent cases, which are useful for BEC based 
interferometers, become much easier with our approach. 

In Chap. 4, I consider a state of iV = l x n bosons that is derived by symmetrizing 
the u-fold tensor product of an arbitrary Z-boson state. The rationale behind this 
approach comes from the BBGKY hierarchy: The errors to the many-particle RDMs 
only weakly affect the few-particle RDMs. Particularly, we are interested in the pure 
state case for l = 2, which we call the Pair-Correlated State (PCS). I show that the 
PCS approach reproduces the number-conserving Bogoliubov approximation when 
depletion is low; moreover, PCS allows one to go to the strong interaction regime 
where the Bogoliubov approximation fails. The normalization factor of the PCS is 
calculated in the large N limit; the few-particle reduced density matrices can then 
be derived by taking derivatives of the normalization factor with respect to the PCS 
parameters. For the two-mode case, these matrix elements are related to the PCS 
parameters by modified Bessel functions. In the basis that the single-particle RDM is 
diagonalized, the two-particle RDMs of PCS has large corrections corresponding to 
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annihilating two particles in one mode and then creating two particles in the other 
mode; although such correlations are crucial for strongly correlated bosonic systems, 
they are always zero in the doublc-Fock ansatz. 

In Chap. 5, the two-site Bose-Hubbard model is used to benchmark the PCS ansatz. 
I find that the error (measured by the trace distance between the two-particle RDMs) 
of PCS to the exact solution is less than two percent over the entire parameter space 
(the error of DFS is roughly 10 times larger). Thus PCS serves as a bridge between the 
superfluid and the Mott insulating phases. More interestingly, for the time-dependent 
case, numerical simulations suggest that the error of PCS does not become larger 
as time increases. I derive both time-dependent and -independent equations for the 
ground state and the time evolution of the PCS ansatz, along with a condition for 
fragmentation. The time complexity of simulating PCS does not depend on N and is 
linear in the Schmidt rank of the two-particle state used to construct PCS. Compared 
to other methods, e.g, the Jastrow wavefunction, the Gutzwiller wavefunction, and 
the multi-configurational time-dependent Hartree method, our approach does not 
require quantum Monte Carlo nor a particularly demanding computational power. 


1.6 Other Work 

In addition to the work on BECs reported in this dissertation, I also took advantage 
of the broad research interests in the Center for Quantum Information and Control 
(CQuIC) at the University of New Mexico to work on a variety of other research topics. 
In the following, I describe briefly some published work that I participated in during 
my PhD study. Typically, my role in these other projects arose from my realizing 
that I could make a contribution, based on my expertise in quantum measurement 
and information theory, to problems that were discussed in the group meeting of my 
supervisor, Professor Caves. 

A quantum linear amplifier makes a small signal larger, so that it can be perceived 
by instruments incapable of resolving the original signal, while sacrificing as little 
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as possible in signal-to-noise ratio. Quantum mechanics limits how well this can 
be done: The noise added by the amplifier, when referred to the input, must be at 
least half a quantum at the operating frequency. This well-known quantum limit 
on deterministic linear amplifiers only constrains the second moments of the added 
noise. In [CCJP12], we derived the quantum constraints on the entire distribution of 
added noise. We showed that any phase-preserving linear amplifier is equivalent to a 
parametric amplifier with a physical state for the ancillary mode; the noise added 
to the amplified field mode is distributed according to the Wigner function of the 
ancilla state. 

A noiseless linear amplifier takes an input coherent state to an amplified coherent 
state, but only works part of the time. In [PJCC13], we bounded the working 
probabilities of probabilistic and approximate noiseless amplifiers and constructed 
theoretical model amplifiers that achieve some of these bounds. Our chief conclusions 
were the following: (i) the working probability of any phase-insensitive noiseless 
amplifier is very small in the phase-plane region where the device works with high 
fidelity; (ii) phase-sensitive noiseless amplifiers that work only on coherent states 
sparsely distributed on a phase-plane circle centered at the origin can have a reasonably 
high working probability. 

Any evolution described by a completely positive trace-preserving linear map 
can be imagined as arising from the interaction of the evolving system with an 
initially uncorrelated ancilla. The interaction is given by a joint unitary operator, 
acting on the system and the ancilla. In [JPC13], we determined the properties 
such a unitary operator must have in order to force the choice of a physical—that is, 
positive—state for the ancilla if the end result is to be a physical—that is, completely 
positive—evolution of the system. Thus Ref. [JPC13] finds the general solution to 
this problem, which reveals a surprising and previously unsuspected structure in the 
joint unitary operators. The problem, in a more restricted setting, arose in our work 
on the general limits to the noise added by deterministic linear amplifiers [CCJP12], 

In quantum optics a pure state is considered classical, relative to the statistics of 
photodetection, if and only if it is a coherent state. A different and newer notion of 


Chapter 1. Introduction 


15 


nonclassicality is based on modal entanglement. One example that relates these two 
notions is the Hong-On-Mandel effect, where modal entanglement is generated by a 
beamsplitter from the nonclassical photon-number state | 1) <g) | 1). This suggests 
that beamsplitters or, more generally, linear-optical networks are mediators of the two 
notions of nonclassicality. In [JLC13], we showed the following: Given a nonclassical 
pure-product-state input to an TV-port linear-optical network, the output is almost 
always mode entangled; the only exception is a product of squeezed states, all with 
the same squeezing strength, input to a network that does not mix the squeezed and 
antisqueezed quadratures. Our work thus gives a necessary and sufficient condition 
for a linear network to generate modal entanglement from pure-product inputs, a 
result that is of immediate relevance to the boson-sampling problem. 

In [Jial4], I derived explicit expressions for the quantum Fisher information and 
the Symmetric Logarithmic Derivative (SLD) of a quantum state that is expressed 
in exponential form p = exp(G'): the SLD is expressed in terms of the generator G. 
Applications include quantum-metrology problems with Gaussian states, including the 
effects of photon losses or other decoherence, and general thermal states. Specifically, 
I gave the SLD for a Gaussian state in two forms, first, in terms of its generator 
and, second, in terms of its moments; the Fisher information was calculated for both 
forms. 

Probabilistic metrology attempts to improve parameter estimation by doing a 
measurement and post-selecting states that are especially sensitive to the parameter in 
question. Such probabilistic protocols thus occasionally report an excellent estimate 
and the rest of the time either guess or do nothing at all. In [CFJC14], we showed 
that such post-selected probabilistic protocols can never improve quantum limits on 
estimation of a single parameter, both on average and asymptotically in number of 
trials, if performance is judged relative to the mean-square estimation error. We 
extended the result by showing that for a finite number of trials, the probability of 
obtaining better estimates using probabilistic metrology, as measured by mean-square 
error, decreases exponentially with the number of trials. 
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Chapter 2 

Basics of Bose-Einstein 
Condensates 


From a certain temperature on, the molecules “condense” without attrac¬ 
tive forces; that is, they accumulate at zero velocity. The theory is pretty, 
but is there some truth in it. 

- Albert Einstein 

Bose-Einstein condensation has attracted many theoretical efforts since its first 
realizations [AEM + 95, DMA + 95]. Perhaps one reason is that BECs are much “cleaner” 
than other many-body systems and thus can be modeled by simple Hamiltonians. 
The “simpleness” of BECs makes possible theoretical investigations of important 
concepts such as nonlinear effects, elementary excitations, and macroscopic quantum 
coherence [P056, Yan62], In this chapter, I review the theoretical frameworks for 
describing BECs, which include the definition of a condensate, the Gross-Pitaevskii 
Equation (GPE) [Gro61, Pit61], and relative phases of BECs [ATM + 97]. In particular, 
the GPE is derived by projecting the evolved many-body state back to the manifold of 
product states; a similar method is also used in Chap. 5 to derive the time-dependent 
equations for the so called Pair-Correlated State (PCS). 1 

1 For situations where depletion (noncondensate fraction) is large, neither the GPE nor 
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In 1924 Einstein received a letter from Bose on the derivation of Planck’s law 
using a new statistics for photons, later known as Bose-Einstein statistics. He soon 
realized the importance of that work, applied Bose’s idea to massive particles, and 
predicted a new phase of matter—Bose-Einstein condensation—where almost all 
the particles condensate to a single quantum state when the temperature drops 
beneath a critical temperature T c . Interestingly, the condensate is a consequence 
of Bose-Einstein statistics, which dramatically suppress the number of low-energy 
excited states, rather than due to attractive particle-particle interactions. Although 
the idea of BEC was initially proposed for noninteracting particles, it also applies 
when there exist interparticle interactions, but great care needs to be taken, because 
the interparticle interactions deplete part of the condensate even at zero temperature. 
To strictly define BECs for interacting bosons at zero or finite temperature, one 
introduces the single-particle Reduced Density Matrix (1RDM), 

P (1) (x|x / ) = <^|rJ; t (x / )iKx)|^>, (2.1) 

where we use upright Greek letters to denote field operators and slanted capital Greek 
letters to denote many-body states. The Penrose-Onsager criterion [P056] states that 
a BEC occurs when the largest eigenvalue of p ^ is of order N, where N is the total 
number of particles. This criterion remains meaningful for interacting systems and 
corresponds to the existence of Off-Diagonal Long Range Order (ODLRO) [Yan62], 
In this chapter, we only consider the case of weakly interacting particles with small 
depletion, in which case the largest eigenvalue of p ^ does approach N. When there 
are interparticle interactions, however, the bosons no longer condense to the ground 
state of the single-particle Hamiltonian, but rather condense to the ground state of 
the time-independent Gross-Pitaevskii equation 

h 0( x ) = (~|^V 2 + ^(x) + 9(N - l)|0(x)| 2 )0(x) , (2.2) 

where 0(x) is the condensate wavefunction, V (x) is the trapping potential, N is the 

number of particles, and g represents the strength of the inter-particle interactions. 

the Bogoliubov approximation works, and one usually needs to rely on methods such as 
quantum Monte Carlo. In Chap. 5, we discuss how to apply the PCS ansatz introduced in 
Chap. 4 to fragmented BECs. Compared to other approaches, PCS is much less numerically 
demanding while capturing interparticle correlations. 
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Since the atoms are very cold, their interaction can be described by the single 
parameter g = Airh 2 a s /rn , where a s is the s-wave scattering length and m is the 
mass per particle. This situation, in stark contrast to the case of liquid helium, 
is a simplification that collisions can be described in terms of the lowest-energy 
scattering length [Fer36, Hua87]. Notice also that in the Gross-Pitaevskii equation, 
the interaction term in Eq. (2.2) is “amplified” IV — 1 times by the Bose statistics. 
Thus, even for weakly interacting dilute Bose gases, the interaction term in the GPE is 
considerable; the effects of interactions are less crucial when temperature is increased 
towards the critical temperature. 


2.1 Derivation of the Time-Dependent GPE 

Generally we need to rely on approximations to deal with interacting many-body 
quantum systems. The inner product between the approximated and the actual 
state, whose absolute value is called the fidelity, is usually considered a measure 
of merit of the approximation. As the dimension of the Hilbert space becomes 
large, however, this inner product goes to zero. One example of this is that the 
inner product between the GP ground state, which is a product state, and the more 
accurate Bogoliubov ground state is quite small. Such a situation urges one to find a 
more useful measure of the merit of approximations that does not go to zero as the 
number of particles becomes large. Often it is the case that the important physical 
properties of many-body systems can be determined by the correlation functions 
(also called Green’s functions) instead of the whole many-body state vector, and an 
approximation may be considered good if it gives the correct low-order correlation 
functions. We emphasis that this criterion is weaker than the inner product criterion; 
the low-order correlation functions of two states can be almost identical, while the 
whole many-body wavefunctions are vastly different. 


Of particular interest are the equal-time 2g-point correlation matrices, i.e., the 
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g-particle Reduced Density Matrices ((/RDM), 

t) = <*'(«) |'t |, (x' 1 )-"l|) , (x',) v|j(x,)- ■ ,ij>( x i) |!?(«)} , 

(2,3) 

where the field operators are time independent (Schrodinger picture). The time 
derivative of the (/RDM is 

ihp (q \t) = (&(t) I ['4 ,t (x 1 ) • • • '4 ,t ( x g) tKxg) • • • ^(xi), H(t)] |^(t)> . (2.4) 

If the Hamiltonian only contains terms up to two-body interactions (quartic in the field 
operators), the commutators [H, rjd(x y )] and [H, 4>(xj)] are at most cubic functions 
of the field operators, which suggests that p ^ is only a function of p^ +1 ^. Generally, 
the fcth time derivative of the (/RDM, d k p^/dt k , is a function of p^ q+k ^ [BBG], and the 
short time evolutions of the low-order RDMs are immune from errors of higher-order 
RDMs. 

Before going further, let us review a general procedure for approximating state 
evolution in a Hilbert space H by a set of ansatz states A C H, where the subset A is 
not necessarily a subspace of H (it can be a submanifold of H). Initially, the state of the 
system is assumed to be in the ansatz subset, | ^(0)) 6 A. After evolving it for a short 
period of time dt, we “project” the evolved state | E(dt)) = exp {—ihidt/K) | 'P'(O)) back 
to A. The “projection” step simply means to find the normalized state | \P A (dt) ) G A 
such that the inner product (\P A (dt) \ & (dt)) is real and its norm is maximized. 
Repeating the same procedure for n = tjdt steps, we find a state | E A (t)) G A within 
A that approximates the actual state | \P(t)) = exp (—iHt/h)\ *^(0)). Note that this 
procedure does not in any way guarantee a large overlap, (& A (t) \ E(t )), between 
the actual state and the approximated state. Nevertheless, at each step, it is the 
best approximation to the evolution within A, and we are satisfied if it gives good 
approximations for the low-order RDMs. The adequacy of the approximation is 
determined by the system Hamiltonian and a judicious choice of the subset A. 

I now discuss how to derive the Gross-Pitaevskii equation using the procedure 
described in the preceding paragraph: moreover, the question why the GP ansatz 
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works so well is investigated. The GP ansatz consists of choosing A to be all of the 
product states; in the second-quantized picture, these product states have the form 


I^sp(*)> = H(t)] l vac ) > 

where aL t \ is the creation operator for the condensate mode, 


4(0 = / ^ T ( x )</ , ( x > t ) rfx = (^|0(^)> = l^ T ) • 


t • 


(2.5) 


( 2 . 6 ) 


Here we introduce a shorthand notation for the integral as a bra-ket inner product 
between a single-particle state and the field operator. The bra-ket notation introduced 
here, though ad hoc , is useful for manipulating the complicated expressions that arise 
as we proceed. We define the error vector of the product ansatz as 


*err(t) ) = 


m 

ih 


Vgoit) ) - I * gP (t) ) 


gp 

N 


(^/T N ~ N d * (t) N 1 ) 1 vac > ’ 


(2.7a) 

VMV ih L “* WJ yi—/> ( 2 - 7b ) 

where d\^ t \ = j r|h(x) 0(x, t) dx = (\Jj | <fi(t )) and the many-body Hamiltonian "H(f) 
takes the form 

U(t)= j ij; t (x)^-^-V 2 + H(x,t))il;(x) + |[-i|; t (x)] 2 il; 2 (x)dx. (2.8) 


We notice that 


H (aj ,) N | vac) = N 4> 


h 2 2 

2^ V + 


V + g(N — 1)|0| 2 0^(aJ,) A ' | vac) 


+ N (N — 1)^| J (4>5 y 4> 2 d^j (4) N 2 I vac) 

- | A(A- !) (ay A I vac) , 


(2.9a) 

(2.9b) 

(2.9c) 


where 4>^(x, t) = ijr (x) — aj,(t)0*(x, t) is the field operator for the noncondensate 
modes, and 


V{t) =9 dx . 


( 2 . 10 ) 


The term Eq. (2.9c) contributes only an overall phase and will be neglected hereafter. 
By letting </>(£) satisfy the Gross-Pitaevskii equation, 

j-2 

ihtj>(x,t) = (-^V J + l"(x,«)+ 9 (JV - l)|0(x,t)| 2 )^(x,t) , 


( 2 . 11 ) 
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we have 


4(t) = (^W)> = ^ 




"ls v2+v,(f)+9(jv “ 1) i 0(i)|: 


m > , (2.i2) 


which cancels the term (2.9a). Thus, the error vector takes the form, 

| j [^l( x ,i)] 2 0 2 ( x ,i)o? x ^ [4 ( thT” 2 | v ac> , ( 2 - 13 ) 

which comes solely from the term (2.9b). This kind of error, which involves two 
particles excited out of the condensate mode, takes the state outside the product 
ansatz (2.5), where only one particle is allowed to be excited. Thus the GP solu¬ 
tion (2.11), by taking into account the term (2.9a), is optimal in minimizing the 
length of the error vector | ^ err (t) )• 


^err (*) ) = 


N(N-l) 
ih y/N\ 


For gN ~ 1, we have 

l^errWI 2 ~ ^ ~ |?4pW| 2 , (2-14) 

which says that the error is not small at all, and the GP ansatz | <^ p (£)) will deviate 
substantially from the actual state. The error of the 1RDM, however, vanishes, 


pS(x, x) = ( y(x') \|>(x) I !?„) + H.c. = 0 , 


(2.15) 


where H.c. means exchanging x and x ; and taking the complex conjugate. 


The error in the two-particle Reduced Density Matrix (2RDM) is 


P (2) ( 

r err \ 


^1? ■ 


Xi,X. 


y 


= < ^gp I ^ t (x'i)^ t (x' 2 ) 4 , (x 2 )4 , (x 1 ) I v e „ ) + H.C. 

= — — (0*(xi)0*(x2) ( vac | aJ r ' 2 rKx 2 )rJ)(x 1 ) | E err ) + H.c.) . 


(2.16a) 

(2.16b) 
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To evaluate Eq. (2.16b), we notice 
ih 

-^= ( vac | a%~ 2 \Kx 2 )\Kx 1 ) | & eiT ) (2.17a) 

= |/ ^-^yy(vac|aJ r_2 iKx 2 )\Kx 1 )[il; 1 L (x)] 2 (ay JV_2 |vac)(ix (2.17b) 

= | / 0 2 (x)(vac|^(x 2 H(x 1 )[^(x)] 2 |vac)dx (2.17c) 

= g J 0 2 (x) ( vac | ^( x i)^±( x ) | vac ) ( vac | 4>(x 2 ) 4> ^ (x) | vac ) dx (2.17d) 

= 9 J 0 2 (x) (<H x i, x) - 0( x i)0*( x )) (<S(x 2 , x) - 0 ( x 2 ) 0 *( x )) dx (2.17e) 

= ^</ , ( x l)0( x 2)(^( x l, x 2) - I0( x i)| 2 - I0( x 2)| 2 +'7/^) • (2.17f) 

Putting Eq. (2.17f) into Eq. (2.16b), we have 

Pe rr (xi,x 2 , Xi,X 2 ) 

= — — (0*(xi)0*( x 2) (vac | aj r_2 il;(x 2 )^(x 1 ) | l^ err ) + H.c.) (2.18a) 

= — ^ ~ 0* (x) )0* (x 2 )0(x 1 )0(x 2 ) 

x (<K x i, x 2 ) - |0( x i)| 2 - |0( x 2 )| 2 + v/g) +H.c. . (2.18b) 

The expression (2.18b) says that fr 2 „ is of order gN 2 (or N) and is negligible compared 
to p ® which is of order N 2 . 

Similarly, for the error in the gRDM, we have 

Peri ( X 1 j ■ ■ ■ > x q ; X 1, ■ ■ ■ , x g) /N(N ~ 1) ■ ■ ■ (N - q + 1) 

= (^*( x i) • • -d*( x g) (vac | a } r ~ f/ cKx f/ ) ■ ■ -^( x i) | ^ err > + H.c.) (2.19a) 

= Jr nt'MM*) E (^( x i’ x fe) _ l^( x i)| 2 “ l<K x fc)| 2 + d/d) + H.c. (2.19b) 

1=1 j ,k=l 


From Eq. (2.19b), we observe that p^l ~ N 1 for q <C vW. This justifies that 
the product ansatz is a good approximation provided that the evolution time is not 
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too long. 2 


2.2 The Relative Phase of Two BECs 

Although it is not possible to attribute an overall phase to a BEC with a fixed number 
of particles, the relative phase of two BECs can be well defined when there is an 
uncertainty in the relative number of bosons. This is linked closely with the well- 
known uncertainty relationship between phase and number [Cav81]. The many-body 
wavefunction is therefore a superposition of states with different numbers of particles 
distributed in the two BECs; such states evolve at different rates due to their different 
chemical potentials caused by the particle-particle interactions. This effect, called 
phase diffusion, degrades the relative phase of the two BECs. In App. A, I discuss 
the equivalence of two seemingly different approaches to phase diffusion. 

Many authors have discussed the relative phase of two BECs from a variety of 
perspectives. For example, Imamoglu [ILY97] and Villain et al. [VLD + 97] showed 
that the phase memories of BECs are lost on a relatively short time scale, in some 
cases vanishing in the large N limit. 3 Preparation and detection of the relative phase 
of two BECs using optical means were discussed in [Jav96, RW97, HMWC98]. Horak 
and Barnett [HB99] showed that coherence between two BECs can be generated 
by just measuring several atoms, but collapsing to the phase state requires order 
N measurements. Zapata et al. [ZSL03] studied the dynamics of the relative phase 
following the connection of two independently formed BECs. Saba et al. [SPS + 05] 
performed the first nondestructive measurement of the relative phase of two spatially 
separated BECs by stimulated light scattering. Chwedehczuk et al. [CPS11] showed 
that the sensitivity of phase estimation by measuring the position of atoms saturates 
the bound set by the quantum Fisher information. 

2 

In practice, the product ansatz probably also works for very long time t, but to prove 
it rigorously is hard. 

3 The coherence time is less than 50 ms for a typical million-atom BEC (with diluteness 

parameter na 3 ~ 10 4 ). 
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Condensed-matter physicist Philip W. Anderson once asked the question: “Do 
two superfluids which have never seen one another possess a relative phase?” This 
intricate question was partially answered by the MIT experiment [ATM + 97], where 
two BECs of sodium atoms are formed separately in a double-well trapping potential. 
The double-well potential is a result of a combination of a magnetic trap and sheet of 
blue-detuned far-off-resonant laser light in the middle. The magnetic trap and the 
laser-light sheet are suddenly switched off, and after about 40 ms time-of-flight, the two 
expanding condensates overlapped and were observed by absorption imagining. Very 
clear fringe patters were observed, and the fringe period was the de Broglie wavelength 
associated with the relative motion of atoms. The centers of the fringe patterns, 
however, were different from shot to shot. This experiment implies that two separate 
BECs possess a relative phase, although this phase might be completely random. 
Many authors discussed this phenomenon using different theoretical approaches; 
see [JY96, NWS + 96, CGNZ96, CD97b, PS08]. Here, I briefly go over the argument 
in [NWS + 96], which makes use of the correlation functions. 


The quantum state of two separated BECs, each with n atoms, can be approxi¬ 
mated by the double-Fock state 


l^dfs) = -7 (o)i) n K) n |vac) , ( 2 . 20 ) 

n! 

where aj) and a R are creation operators of the left and the right single-particle states 
'0 L (x) and One simple way to explain the experimental results is by calculating 

the correlation of particle densities [NWS + 96], 


q(xi,x 2 ) 



vac 


«. OdH 1 ’ { x 2 ) (°r) n ( 4 ) n 


vac 


= n(n - 1) \A{xi)\ 2 \^l{x 2 )\ 2 + n(n - 1) |^ R (x 1 )|-|^ a (a; 2 )|“ 

+ n I^lOdO^rOd) +^l(^i)^ r (^ 2 )| 2 

~ n 2 ( I^lOd)! 2 + IV'r^i)! 2 ) (I^ l (^ 2 )I 2 + I^r(^ 2 )| 2 ) 

+ n + c.c.) , 


4 


In this section, we assume the BECs are quasi one-dimensional. 


(2.21a) 

(2.21b) 


(2.21c) 
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where the Fock terms in Eq. (2.21c) describe the interference of the two BECs. With 
the correlation function, we have the expectation of the squared norm of the Fourier 
component of the particle density distribution, 


e lkx 4>t (x)r[)(a;) dx 


= I {^\x 1 )^{x^\x 2 )^(x 2 )) dx x dx 2 (2.22a) 


= f e ik{xi ~ X2) (^ f ( 

= 2n + [ e lk(<Xl ~ X2> g(x 1 ,x 2 ) dx x dx 2 , 


(2.22b) 


where the constant 2 n comes from different ordering of the creation and annihilation 
operators. Neglecting this constant, we define 


Qk 


e lk ^ Xl X2> g(x 1 , x 2 ) dxidx 2 . 


(2.23) 


Consider the case where the left and right states are plane waves in an interval of 
length L, i.e., = e lk ° x /\[L and = e~ lk ° x /\/L- 1 we have the following for k ^ 0, 

2 r- 



e tk ^ Xl X2) (e 2lko(Xl X2) + c.c.) dx ± dx 2 = n 2 (6 kj2ko + d k) _ 2ko ) 


(2.24) 


The quantity g k is nonzero for k = ±2 k 0 , and thus one only hnds a fringe pattern with 
period i r//c 0 ; the fluctuation of the contrast of the fringe patterns is zero from shot to 
shot, a result that can be proved by calculating higher-order correlation matrices. 
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Chapter 3 

Number-Conserving Bogoliubov 
Approximation Using Extended 
Catalytic States 


With every passing year, BEC proves that it still has surprises left for us. 

- Eric A. Cornell 

In this chapter, we consider the ground state and dynamics of a dilute-gas BEC of 
N bosonic atoms trapped in an arbitrary external potential. In order to describe 
how interparticle correlations modify the Gross-Pitaevskii equation, we go to the 
next level of approximation, the Bogoliubov approximation. The Bogoliubov approxi¬ 
mation [Bog47, Fet72, Hua87] is important for several reasons: (i) it tells when the 
Gross-Pitaevskii (mean-held) approach begins to break down; (ii) it describes small 
deviations from the Gross-Pitaevskii equation and can be used to study the stability 
of a BEC; (iii) it enables the calculation of how impurities change the behavior of a 
BEC; and (iv) it is useful for studying phase coherence between BECs. 

Conventionally, in the Bogoliubov approximation, the condensate is treated as 
a small perturbation of the state where all the bosons occupy a coherent state of 
a particular condensate mode (i.e., a single-particle state). When particle loss is 
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negligible, however, the real condensate is much closer to a number state than to a 
coherent state (see Fig. 3.1). Since a coherent states has a well-defined phase, the 



Figure 3.1: Phase-space representations for (a) a coherent state with complex 
amplitude a and (b) a Fock (number) state with particle number N. A Fock 
state is distributed phase-symmetrically on the phase space, and no definite phase 
can be attributed to it; in contrast, a coherent state has a well defined phase. 

conventional Bogoliubov approximation breaks the f/(l) symmetry possessed by the 
condensate; consequently, a fictitious Goldstone mode [Nam60, Gol61] is present in 
the Bogoliubov Hamiltonian (see Fig. 3.2). Because there is no restoring force on 
the Goldstone mode, the Bogoliubov ground state is not well defined; worse, the 
Goldstone mode causes the condensate state to deviate linearly in time from a single 
condensate in a coherent state (i.e., this is a secular deviation, not an oscillation). This 
problem is particularly pesky when the condensate is in a trapping potential, where 
the Goldstone mode is a mixture of the condensate mode and modes orthogonal to it 
and thus cannot be removed easily. The way to get rid of the unphysical Goldstone 
mode is to adhere to the fact that the condensate has a fixed number of particles, 
i.e., by using a Bogoliubov approximation where particle number is conserved. 

Many authors have already considered the number-conserving Bogoliubov approx¬ 
imation. Girardeau and Arnowitt [GA59] were the first to propose a theory for the 
ground state and excited states of many bosons based on a particle-number-conserving 
(iV-conserving) formulation of the Bogoliubov quasiparticles and quasiparticle Hamil¬ 
tonian. Gardiner [Gar97] introduced a somewhat similar approach to Girardeau 
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Figure 3.2: The breaking of U(l) symmetry of the condensate wavefunction causes 
a Goldstone boson in the Bogoliubov Hamiltonian. Because there is no restoring 
force on the Goldstone mode, the Bogoliubov ground state is not well defined. 

and Arnowitt’s, but with an emphasis on the time-dependent case; Gardiner et 
al. [GZBD97] then applied this approach to the kinetics of a BEG in a trap. Castin 
and Dum [CD97a, CD98, CasOl] gave a modified form of the Bogoliubov Hamilto¬ 
nian where the terms that break the U( 1) symmetry are removed by a projection 
operator. Sprensen [Spr02] generalized the Castin-Dum result to the two-component 
case. Gardiner and Morgan [GM07] considered an expansion in powers of the ratio 
of noncondensate to condensate particle numbers, rather than in the inverse powers 
of the total number of particles. Leggett [LegOl, Leg03] used a number-conserving 
BCS-like ansatz to discuss the properties of the ground state of a homogeneous BEC. 
The so-called truncated Wigner approximation [SOP + 98, SLC02] is an equivalent 
way to implement the number-conserving Bogoliubov approximation in a phase space. 

A number-conserving Bogoliubov approximation yields qualitatively different 
results from one that fails to conserve particle number; among these are the following. 
Villain et al. [VLD + 97] showed that the collapse time of the phase of a BEC is 
relatively short and, in some cases, vanishes in the limit of a large number of atoms. 
Danshita et al. [DEYK05] investigated collective excitations of BECs in a box-shaped 
double-well trap. Trimborn et al. [TWK08, TWK09] demonstrated that artificial 
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number fluctuations lead to ambiguities and large deviations within the Bose-Hubbard 
model. Oles et al. [OS08] predicted considerable density fluctuations in finite systems 
close to the phase-separation regime. Schachenmayer et al. [SDZ11] studied the 
collapse and revival of interference patterns in the momentum distribution of atoms 
in optical lattices. Gaul and Muller [GM11] studied BECs in spatially correlated 
disorder potentials. Billam et al. [BMG13] derived equations of motion describing 
the coupled dynamics of the condensate and noncondensate fractions. 

We return to the number-conserving Bogoliubov approximation here to provide a 
particularly transparent method of deriving the relevant equations. Our approach 
to a number-conserving Bogoliubov approximation is to “encode” the many-body 
wavefunction of the BEC in the iV-particle sector of a state we call an Extended 
Catalytic State (ECS), by which we mean a coherent state for the condensate mode 
and a state to be determined by the dynamics for the orthogonal modes of the atoms. 
Using a time-dependent interaction picture, we move the coherent state to the vacuum, 
thus making all the field operators formally small compared to TV 1 / 2 . The resulting 
Hamiltonian can then be organized by powers of N^ 1 ^ 2 . Requiring the terms of order 
N 1 ^ 2 to vanish in the interaction-picture evolution equation gives the GP equation 
for the condensate wavefunction. Going to the next order in the evolution equation, 
N°, we derive equations equivalent to those found by Castin and Dum [CD98] for a 
number-conserving Bogoliubov approximation. In contrast to other approaches, ours 
allows one to calculate the state evolution in the Schrodinger picture, and it also has 
advantages in considering higher-order corrections and extensions to multi-component 
cases. 


3.1 Encoding the State of a BEC 

What is true in quantum optics, that coherent states are easier to deal with than 
number states, is often true elsewhere. Indeed, the usual mean-field approximation 
to BEC evolution is based on the assumption that the BEC is in a coherent state 
of a condensate mode [PS03]. A problem with this approach is that the number 
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of particles in a BEC is usually fixed, whereas coherent states are superpositions 
of states with different numbers of particles. A related problem is that assigning 
a coherent state to a BEC breaks its phase symmetry, thus causing problems in 
developing the Bogoliubov approximation. 

Our philosophy for dealing with a BEC that has N particles is to extend the BEC 
state | \P N ) to a state | \P ecs ), for which the condensate mode is in a coherent state, but 
the IV-particle sector is the same as ) within a constant. Consider an arbitrary 
state | W N ) with N particles, for which we have the number-state decomposition, 

N 

\E n ) = Y j \N-M),®\C m ) ± , (3.1) 

M =0 

A/j_ | C M )j_ — M | i? M )j_ , (3.2) 


where the kets labeled by 0 and _!_ apply to the condensate mode and to all the modes 
orthogonal to the condensate mode, respectively. The operator A f± is the particle- 
number operator for the orthogonal modes. The state | ) for the orthogonal 

modes, which has M particles in the orthogonal modes, is not necessarily normalized. 
The key to our approach is that the state (3.1) can be written as 

N 


v N ) = e '«' a/2 V klAzih! Vn ( | a )o »| n M h 


M =0 


a 


\a\ 2 /2 VN~ ] - sp, (| \ ^ ( M / {N ~ M)\ x ^ 

— e — n~ I I a )o ® ( 2.^ a Y —An— I ^ M ) 

a \ ^ m= o ' ' 


= e H 2 /2 


a 


N 


Vn{\ oi )o <8 | Cl) 


(3.3a) 

(3.3b) 

(3.3c) 


where Vn is the projection operator onto the IV-particle sector and 

N 


M =0 




N\ 


(3.4) 


is an (unnormalized) state of the modes orthogonal to the condensate mode. 
We now introduce the extended catalytic state , 


S'ec) = I a)o ® I a)± . 


(3.5) 
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which is related to the physical state by 

I &N ) — e |a|2/2 — V N \ tl/ecs ) . (3.6) 

a 

The extended catalytic state is a direct product of a coherent state |a) 0 in the 
condensate mode and an unnormalized state | fi)±_ of the orthogonal modes. Notice 
that once a is specified, the extended catalytic state has a one-to-one correspondence 
with the physical state. The structure of the extended catalytic state allows us to 
study the dynamics of a BEC in the Schrddinger picture, and we will see that the 
structure is preserved throughout the evolution in the Bogliubov approximation. 

For a pure condensate with no depletion of the condensate mode, the modes 
orthogonal to the condensate mode are in vacuum, and the overall state has the form 

| \P N ) = | N) 0 ® | vac) ± . (3.7) 

In this case we have 

l^ecs) = |a)o® |vac)_L . (3.8) 

Generally one expects that a dilute-gas BEC has a state close to that of a pure 
condensate, in which case the noncondensate state | fl )j_ is close to the vacuum; 
we want to develop an approximate description based on this expectation. To do 
so, notice that the encoding into an extended catalytic state works for any value 
of a. In other words, one has the freedom to choose a at will; after the projection, 
all values of a yield the same physical state. Nonetheless, we stick to the choice 
| or | = iV 1 / 2 , for the reason that we make approximations in deriving the dynamics 
of | E ecs ) and the projection into the iV-particle sector can amplify the errors due 
to these approximations. To keep these errors under control, we need the number 
distribution of the coherent state to be centered at the actual atomic number N. The 
phase of a is yet another matter, which we discuss further below. 

The BEC Hamiltonian conserves particle number and thus commutes with the 
particle-number operator. As a consequence, the evolution operator U(t) commutes 
with Vn, allowing us to move the evolution operator through the projection into the 
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IV-particle sector so that it acts directly on the extended catalytic state: 


V N {t))=U{t) \V N (0)) (3.9a) 

= e \°\ 2 /2^M U (t)V N |^ ecs (0)) (3.9b) 

a 

= eW 2/2 ^V N U(t)\V ecs (0)) (3-9c) 

a 

= e'“ |2/2 V N \ *„{*)) • (3.9d) 


To find | !^v(t)), one solves for | \P ecs (t) ) and then projects into the iV-particle sector. 

The first step in developing our approximation is to go to an interaction picture 
in which the condensate mode is displaced from a coherent state to vacuum. To do 
this, we start with a time-dependent condensate mode defined by the single-particle 
state | (f)(t )), which has wave function 

0(x,t) = (x|0(t)) . (3.10) 

The annihilation operator for the condensate mode is related to the Schrodinger- 
picture field operator i|>(x) by 

= J 0*(x,f)rKx)dx = (<f>(t) |\p) = (V !</>*(£)) • (3.11) 

Here, in the final two equalities, we introduce a shorthand notation for the integral 
as a bra-ket inner product between a single-particle state and the field operator. The 
creation operator for the condensate mode is 

4(0 = f 4(x)0(x,f)dx = (tp|0(t)) = |rl4 . (3.12) 

The bra-ket notation introduced here, though ad hoc , is useful for manipulating 
the complicated expressions that arise as we proceed, more so once we get to the 
two-component case in the next section. The annihilation and creation operators 
have two different bra-ket forms, both of which we use in our treatment. The field 
operator can be written as 


<l>(x) = a m </y(x, t) + t|)j.(x, t) , 


(3.13) 
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where 4> ± (x, i), the held operator with the condensate mode excluded, can be expanded 
in terms of single-particle states orthogonal to j 0(t)). In the Schrodinger picture, 
i0_i_(x, t) has a time dependence because the condensate mode is changing in time. 
Notice that in terms of our shorthand notation, we can write 

I MM = |M>} - \ (/>}a<t> = |M>) - |0)(0|M>) , (3.14a) 

('4’j.I = (4 , l — a 0<0l = (4 , l — ('M , l0)<0l • (3.14b) 

The extended catalytic state for a pure condensate in the time-dependent conden¬ 
sate mode | 0(f)) is 

X>(a, 0(f)) | vac) = | a, 0(f) ) 0 <8> | vac) ± , (3.15) 

where the displacement operator Z>(a,0(f)) for the condensate mode, which we 
usually abbreviate as T>(t), is dehned as 

V(a, 0(f)) = V(t) = exp [aa\ {t) - a*a^) , (3.16) 

The state (3.15), which describes a pure condenstate with no depletion, is the one we 
perturb about in developing our approximate description. 

We can now introduce the desired interaction picture as the one where the 
condensate mode is displaced to vacuum; i.e., states transform to 

_MM_ 

I MM) = V'(a,<j>(t))\* ecs (t)) =P t (a,0(t))W(t)P(a,0(O)) | Mt(0) > , (3.17) 

where 7/ int (f) is the evolution operator in the interaction picture. The Schrodinger- 
picture evolution operator U (f) obeys the Schrodinger equation 

ihU(t) = H(t)U(t) , (3.18) 

where EL{t) is the (possibly time-dependent) BEC Hamiltonian. The time depen¬ 
dence of the condensate wave function, which enters into the displacement operator 
27(a, 0(t)) through the annihilation and creation operators, and aM, is to be 
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determined. The evolution operator U mt (t) obeys the equation 


ihu iDt (t) = ihv\t)u(t)v(o) + ihv ] (t)u(t)v( o) 

= (ih&{t)v{t) + v\t)n{t)v{t)}uM{t) ■ 

The time derivative of the displacement operator is 

t>\t) = y 

= 4 |a|2 j (V* a * 

= — a e a ^ a\ e~ a “ 0 j 1^(t) 

= (a* a ,$ - a dj, - l«| 2 [<ty, 4]) (*) 

= — a 4 — |a|“ (0 | 0)j T>\t) . 

Putting the above expression into Eq. (3.19b), we have 

ihU int (t) = Hintit) U mt (t) , 


(3.19a) 

(3.19b) 

(3.20a) 

(3.20b) 

(3.20c) 

(3.20d) 

(3.20e) 

(3.21) 


where the interaction picture Hamiltonian reads 

■Hint(t) = ~ih ^ \a\ 2 ( 4>(t) | ) +«4(0 ~ a *a<Kt)) +^{t)H{t)V(t) . (3.22) 

Equivalently, we have 

^ I #int(f) ) = ^int(^) I ^intO) ) ■ (3-23) 


In the interaction picture the field operator takes the form 

T>\t) 4>(x) V(t) = i[>(x) + a 0(x, t) . (3.24) 

An expansion of % int (£) in powers of l/|a| = 1/N 1 ^ 2 is a good approximation as long 
as the field operator is small relative to the interaction picture, i.e., more formally, as 
long as the one-particle density matrix is small in the sense that 

Ant(x, x) = ( W int | rjd(x') 4>(x) | W int } ~ N° . 


(3.25) 
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This requirement is satisfied as long as the system is a condensate. In second-quantized 
form, the model Hamiltonian for the BEC is 

f (V(x)(-^v 2 + y(x,t))\|>(x) + |rjj T (x)ii) t (x)rKx)rKx)^ dx , (3.26) 

where the first term is the second-quantized Hamiltonian for particles trapped in a 
potential V(x, t) and the second term represents the two-body scattering energy. The 
only explicit time dependence in the Hamiltonian (3.26) comes from a possible time 
dependence in the trapping potential V(x,t)- For our expansion in powers of l/|a| 
to work, we must have that g\a\ 2 is of order iV°; i.e., g is of order N _1 . 


Going to the interaction picture, we have 




= -ih ^ \a\ (0(t) | </>(£)) + aa\ {t) - a* d^ t) J + V\t) H(t) V(t) 

■-2 


~ a 


■ p 9 h 2 , T/ J 9 I 121 / 1 2 

- ,h Wt~ 2^ V +l/+ 2 |o l 1^1 


dx 


+ J i]; t — ^-V 2 + ^ + 5 , |q;| 2 |0| 2 ^ (f> dx + H.c. 


(3.27a) 

(3.27b) 

(3.27c) 


V ( " + ^ + + f 0 2 + («*) 2 ^ (0*) 2 ) 


dx , 
(3.27d) 


where we neglect terms of order N~ 1 ^ 2 or less. The c-number term (3.27b), of order 
N, is, in the time-independent case, the mean-field energy of the BEC. By requiring 
the linear term (3.27c), of order N 1 ^ 2 , to vanish, we get 

ih<p{x, t) = V 2 + V(x,t) + g\a\ 2 \(j)(x,t)\j </>(x,t) , (3.28) 

which is the celebrated Gross-Pitaevskii equation. Thus the structure of our approach 
is clear: By going to the interaction picture, the mean-field evolution is removed, and 
then by neglecting the terms of higher order than N°, we are left with the quadratic 
Bogoliubov Hamiltonian 


qj = 
ti'bog 


V(-|^v 2 + h + 2sM 2 M 2 )^ 

+ |(a 2 r);W + (a*) 2 xbrK^) 2 ) 


dx , 


(3.29) 
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To reveal the symplectic structure of the Bogoliubov Hamiltonian, it is useful to 
write it in the matrix form 



^Hgp +g\a\ 2 \(/>\ 2 ga 2 4> 2 \ 

v g{(y*f{(t>*f + g\a\ 2 \4)\ 2 J \|V)/ 


(3.30) 


where the colons denote normal ordering of annihilation and creation operators and 
the GP single-particle Hamiltonian takes the form 

tfgp it) = -|)v 2 + Vit) + g\a\ 2 \<f,(t)\ 2 . (3.31) 

Notice that the normal ordering has an effect only on the lower-right corner of the 
matrix. 


Like all the matrices of symplectic structures in this dissertation, we denote the 
2x2 matrix in Eq. (3.30) by a special sans serif character: 


H 


bog 


^ H gp + g\a\ 2 \(j)\ 2 

2 i2 \ 

ga 0 

y S(COV ) 2 

H gP + g\&\ 2 \(t>\ 2 j 


(3.32) 


As shown by Lewenstein and You [LY96], H hog has a nilpotent subspace, where 
phase diffusion takes place. Such phase diffusion is not physical, but rather is a 
consequence of the arbitrary phase assigned to the condensate wavefunction, i.e., 
to a. This problem was addressed by introducing the so-called number-conserving 
approaches [GA59, Gar97, CD98]. Particularly in the work of Cast-in and Dum, a 
systematic expansion of the field operators was used in deriving the equations for the 
number-conserving approach. The aim is to eliminate the artificial nilpotent subspace 
that gives rise to the phase diffusion. Here we solve the same problem by introducing 
an additional contribution to the Hamiltonian, an auxiliary Hamiltonian T{t\ which 
does not affect the Y-particle sector of | \P ecs (t )) and thus keeps the physical state 
| & N (t)) unchanged, i.e., 


VNHt) I ^ecs (t) ) — 0 . 


(3.33) 


With this term tF(t), we can solve the phase diffusion problem by eliminating the 
nilpotent subspace of H hog . 
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To determine the form of J-'(t), we must go to the Bogoliubov level of approxima¬ 
tion, but for now let us suppose 7F(t) takes the form 

Ht) = -^(Af-N ) 2 + (aa\ it) +Af ± (t)-N) T ± (t) + (a*a m -N) P L (t) . (3.34) 
Here 

AT = j ij/(x)4>(x) dx (3.35) 

is the total particle-number operator, and A f± = Af — is the particle-number 

operator for all the modes orthogonal to the condensate mode, i.e., the depletion 
number operator. The time-dependent parameter rj(t), which is to be determined, is 
of order iV _1 . The operator J r ± , also to be determined, is of the order iV~ 1//2 , is a 
linear function of the annihilation and creation operators of the modes orthogonal to 
the condensate mode, and thus commutes with a^ and aj,. The first term in Eq. (3.34) 
clearly satisfies Eq. (3.33). For the other two terms, we have 

"P tv {ota}^ + Af± — iV) J~_\_ | oi )q ® | C )j_ = — iV) | ot )o ® ±\ f?)j_ — 0 , (3.36) 

(< y*a 4 - N)P l | a ) 0 ® \ Q ) ± = (|a| 2 - N) \ a ) 0 ® P ± \Q ) ± = 0 , (3.37) 

where in the first equation we use aaj, + Af± — N = a^(a — a^) + AT — N. As long as 
the condensate mode stays in the coherent state with amplitude a , these two terms 
do not affect the physical state | ^(t )). We show that the condensate mode does 
remain in a coherent state at Bogoliubov order at the end of this section. 

An astute reader might object at this point that the auxiliary Hamiltonian (3.34) 
is not Hermitian. This is not a problem at Bogoliubov order, however, because 
the only nonHermitian term in 7F(t) is Af±J r ±, which, being of order iV -1//2 , can be 
neglected in the Bogoliubov approximation (order N°). Going now to the interaction 
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picture, we have 


J- int (t)=I ^(a,<j>(t))T{t)V(a,<l>(t)) 

(3.38a) 

n 

2 

( J [a*(p* + rjj' 1 ’) (a4> + "4>) dx — N^j 

(3.38b) 

+ 

(cr + o^) + A/j_ — N^j J~ j_ + (a + n^) — IV ^ j) 


= --( 
2 ' 

cu20 T cr ci0 T A/”^ T (ttoj), T A/j_) T cr o^ 

(3.38c) 

V 

_ 

2 

^ |Of | ((20(20 + 0-0(20) + tt Q'0O0 + (o; ) O0O0^ + 0(00 + Of 

CLfj, j 



(3.38(1) 


where the identity |a|“ = N is used to cancel several terms and where terms of order 
N ~ 1 ^' 2 or smaller are neglected in the final form (recall that interaction-picture held 
operators are order N°). Adopting this final, approximate form, we have 




int 


( t ) — —— ^ 2\a\ + O0O0 + (<a ) cl^J -[_ — — |o;| 


(3.39) 


where we normally order the creation and annihilation operators of the condensate 
mode in preparation for incorporating J r int into the main Bogoliubov Hamiltonian. 
This normal ordering introduces a c-number term, — ?; |cr| 2 /2, which could be important 
as a second-order correction to the condensate energy, but which only adds an overall 
phase to the evolving quantum state. Thus we neglect this term henceforth. The 
modified (number-conserving) Bogoliubov Hamiltonian then takes the form 


n 


neb 


n 

n 


bog •' int 

_n( 

b °g 2 ^ 


ol |2 t , 2 t t . / *\ 

2 \a\ o, cl^cl^ j 



(3.40a) 

T QtQj]p J~ l T Ol (20 • 

(3.40b) 


To eliminate the phase diffusion, we choose 
’nit) = g J\(j)(x,t)\ A dx (3.41a) 

= g{m\ 2 \^>) = g(<f>*m 2 \d> m )= 0<**i(^)v>=0<^v> (3.4ib) 
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tF ± (t) = tF^t) =—ga* J 0 *(x, £)| 0 (x, t)| 2 i|)_ L (x, t) dx + H.c. ( 3 . 42 a) 

= rj(t)a *— ga* J 0 *(x, i)| 0 (x, f)| 2 ij>(x) dx + H.c. ( 3 . 42 b) 

= ^«*(( 0 || 0 | J | 0 )( 0 |^) - ( 0 || 0 h’|^)) + H.c. , ( 3 . 42 c) 

where i|)j_(x, t) is the held operator of Eq. ( 3 . 13 ), which has the condensate mode 
excluded. It is now a tedious calculation to show that 

77 d f r> 12 t . 2 t t | / *\2 \ 

•Mnt — ~ a <t> a <t> + a <t> a <f> + ( a ) a (j> a <f> J 

/ x /■ \ ( 3 . 43 a) 

- (|a | 2 4 + («*)%) j 0 *(x, t)| 0 (x,f)| 2 \Jj(x) dx + H.c.J 

= f(|a| 2 <l|)| 0 )( 0 |W 2 |'#'>( 0 l't>) 

+ (a) 2 {ip* | f > {f | (f) 2 | ^ 14> > + H.c.) 

~ s (\ a \ 2 ('i , \’t’)(<t>\\<t>\ 2 \'i’) + («*) 2 <l|> t l(<i*) 2 | 4 > > + H.c.) . 

( 3 . 43 b) 


Translating this into matrix notation, we get 





where the symplectic matrix is 


F int = 9 


1 l«l 2 (QI0l 2 Q - H 2 ) 
^(a*) 2 (Q*(^) 2 Q-(0*) 2 ) 


a\Q<j?Q'-<j?) N 

«I 2 (<?*I0I 2 Q*- W 2 ) / 


Here 


(3.44) 


(3.45) 


Q(t) = 1 - P{t) = t - I 0(f) )(0(t) I (3.46) 

is the projector onto the single-particle space orthogonal to the condensate mode, 
with Q*(t) = 1 — | 0*(f))(0*(t) \. In Eqs. (3.41b), (3.42c), and (3.43), we use the 
bra-ket notation, which is the easiest way to carry out the algebraic manipulations 
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and which also generalizes straightforwardly to the two-component case considered in 
the next section. The modified (number-conserving) Bogoliubov Hamiltonian (3.40a) 
now reads 


'Hncb=2-((' l H (ijhlj^ncb 


*1>) 


(3.47) 


with 


/ 


Hneb = H 


bog 


Fint ~ 


H gp + g\a\ 2 Q\(j)\ 2 Q 


\ 9 (a) 2 Q*(<f>*) 2 Q H gp + g\a\ z Q*\^Q*J 


got 2 Q4> 2 Q* ^ 


| / 12 / 0 * 


(3.48) 


We return now to showing that the condensate mode remains in a coherent state 
under the evolution at Bogoliubov order. This is a bit different from saying that 
the condensate mode is decoupled from the orthogonal modes, as would be the 
case if we were considering the time-independent ground state of the condensate. 
Rather, because | </>(t)) is changing in time, what we need is that the coupling to the 
orthogonal modes is of just right sort to maintain the condensate mode in a coherent 
state. 

To show this formally is quite easy. We begin by noting that the interaction-picture 
evolution equation at Bogoliubov order, 

I ^int(t) ) = Uncb(t) I #int(*) ) » ( 3 - 49 ) 

allows us to write 

ih ^ I ^ int )) =inil 4> I ) + a^neb I l// int ) (3.50a) 

[iha^ [77 nc b, Q><f>\ ^ | ^int } T ^ncb(®c/>| ^int )) ■ (3.50b) 

If we could show that the first term in Eq. (3.50b) vanished, then we could conclude 
that if a 0 ( O ) | 7/ mt (0)) = 0 , i.e., the condensate mode starts in vacuum in the interaction 
picture, then a^ t] | ^ int (i)) = 0 , i.e., the condensate mode remains in vacuum at 
all times in the interaction picture. That the first term in Eq. (3.50b) vanishes is 
thus what we mean by saying that the time-dependent coupling to the orthogonal 
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modes is of the right sort. Showing this involves using the GP equation (3.28), the 
commutator [i|^(x), a^] = —0*(x), and the fact that the only term in Eq. (3.47) 
that contains the annihilation and creation operators of the condensate mode is 
(if) \H gp (t) ip ) = : (r{V \H gp (t)\ rjk ): . Pulling all this together, we get 

= ih(<j)(t) |tp) = \H gp (t)\ ij>) 

= [(tP l^gpWI^), = [H nch (t), a m \ , (3.51) 

which shows that the first term in Eq. (3.50b) vanishes. Equivalent to the statement 
that in the interaction picture the condensate mode stays in vacuum if it begins in 
vacuum is the statement that in the Schrodinger picture the condensate mode is 
always in the coherent state T)(a, 4>{t )) | vac) = | a, (p(t )). As a result, Eq. (3.33) is 
always satisfied, and the auxiliary Hamiltonian T does not affect the physical state. 

Notice that this means that in the extended catalytic state, the condensate mode 
is never entangled with the other modes. When we project the extended catalytic 
state to the IV-particle sector to obtain the physical state of the BEC, however, 
entanglement makes its appearance. 

We can see the effect of the Bogoliubov Hamiltonian (3.47) on the condensate 
mode more directly by dividing the field operators in the Bogoliubov Hamiltonian into 
the contribution from the condensate mode and the contribution from the orthogonal 
modes, as in Eq. (3.13). The condensate mode does not notice the terms with 
projection operators Q and Q* in them, so we get 

■Hncb = 4^(01^10) + a\((j)\H gp \^ ± ) + (x|>_L \H gp \(j))a ( j ) + n nch± (3.52a) 
= a\,((j)\H gp \^) + (tpj. | H gp \<f>)a ( / ) + % ncb ± > (3.52b) 


where 


n 


ncb_L 




(3.53) 


is the Hamiltonian for the orthogonal modes, in which we can use H hog instead of 
H nch because the projectors Q and Q* have no effect. Using the GP equation (3.28), 
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we can rewrite Eq. (3.52b) as 

hincb = -iha\a^ + ih(i\) ± \<j))a (j) + 'H ncb_L ? (3.54) 

from which we can immediately verify the commutator identity (3.51). 


3.2 The Two-Component Case 

In the preceding section we discussed how to derive the number-conserving Bogoliubov 
approximation for a single-component BEC by going to an interaction picture where 
the condensate mode is displaced to vacuum. In this section we show that it is a simple 
task to generalize our method to multi-component BECs. We do the two-component 
case as an example, but the generalization to many components is straightforward. 
In the two-component case the condensate wavefunction, which is generally a single- 
particle state that is entangled between the translational and internal degrees of 
freedom, takes the form 

I 0(*)) = - I 0<r(f)) ® I 0-) (3.55a) 

a 

= i (aq(t) | (t) ) (8) | 1) + a 2 (t) \ 0 2 (t)) 8 | 2) ) , (3.55b) 

where a labels the hyperfine levels, taking on values 1 and 2 in the two-component 

case, and where |aq| 2 + |a 2 |" = l a | 2 = N, with N being the total number of particles. 

The states 11) and | 2) are internal states of the bosonic atoms, which we refer to as 
hyperfine levels because that would be a typical situation. In the subspace spanned 
by | 0!) 8) | 1) = | 1, 0i) and | 0 2 ) 8 | 2 ) = | 2, 0 2 ), the single-particle state that is 
orthogonal to the condensate mode is 

I 0(0)) = \ ( a* 2 (t) | 0i(t)) «) 11) - al(t) | 0 2 (t)) 8 | 2 ) ) . (3.56) 

Notice that 

= I'M*)}® 11) = (^I^W> + —1#)>) . (3.57a) 

\ a a J 

= \Mt))®\ 2 ) = (^\m)~—\m)) • (3.57b) 

\ a a J 
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The field operator that destroys a particle in internal level o at position x is 
t|) ct (x). In our shorthand bra-ket notation for field operators, we have 

*M X ) = ( x I MV ) = ( x 1) • (3-58) 

In the final form we extend our notation by introducing a total field operator 

|M>) = > (3-59) 

cr 

which is a spinor field operator, including both spatial and internal degrees of freedom. 
It gives the hyperfine-level field operators according to (cr 1 1 |)) = | r|v); notice that 
since 4>£.(x) = (t|v| x ) — (xp | cr, x ), we also have (ip | cr) = (|. The spinor 
representation is 

■0( x ) = ( X I0) = X^( x )l a ) ' ( 3 -60) 


The annihilation and creations operators that destroy or create a particle in 
internal level o and spatial wave function 0(x) are 



,i\) J 

^V( x 

o*M x 

)dx = (0 |rJO 

= 

,01^), 

(3.61a) 


a U = J 

^0( x ) 

^t( x ) 

d*=(-l\>a 10) 

= 

k,0) • 

(3.61b) 

The 

annihilation operators 

for the states 

0) and 0) are thus 


Clff) 

= (0 1 tP ) 

1 

* 

a 


L | Tpl ) + a *2{ 02 

1^2)] 

1 — + a 2 a 2,</> 2 ) ! 

(3.62) 

CL(f> 

= (0 1 tP ) 

4 ( 

a 2 ( 0i 

1 "*Pl ) — «l(02 1 

^ 2 )) 

= — (a 2 — a l a 2,<t> 2 ) i 

(3.63) 


where is the annihilation operator that destroys a particle in hyperfine state o 

with spatial wavefunction 0 CT (x, t). The field operator for the atoms in hyperfine 
level cr can be written as 


^K( x ) = + , 


(3.64) 
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where i|v±( x , t) excludes the mode with wave function 0 cr (x, f). The total field 


operator can be written in a variety of forms, 

\A>) = ® lO + l'M*)) (3.65a) 

cr 

= Hh)\ 0(f )) + | ) (3.65b) 

= a^t) I </>(*)) + I (t )) (3.65c) 

where 

\^(t)) = ^±(t))W) (3.66a) 

(7 

= I A> ) - I 4>{t)) ~ a#t) I 4>{t)) (3.66b) 

= |tp) - (I <j>(t)){(/>(t) | + | 0(f))(0(f) |)|^) (3.66c) 


is the total field operator with modes | 0) and j 0) removed and 

I'M*)) = s *(t)l0(f)> + l'M*)> = I'l , }-a*(t)l0(*)) = I ^)-|0(*) )(</’(*) I ) (3-67) 

is the total field operator with only the condensate mode removed. By using our 
bra-ket shorthand, all the manipulations for two components can be made identical 
to that for a single component. 

Just as in the single-component case, we perturb about the extended catalytic 
state for a pure condensate that is in a coherent state for the condensate mode: 

T>{a , 0(f)) | vac) = | a, 0(f) ) 0 ® | vac )j_ , (3.68) 

The physical state is obtained by projecting into the iV-particle sector. 


In the two-component case the model Hamiltonian for the N atoms is 


m) = S/’h(-^ v2 +’ / ^)K‘ ix+ E 'MM* 

(7 G <T,T 

dx . 


(3.69) 


The diagonal terms of the Hermitian matrix ficu CTT give the energies of the internal 
levels, and the off-diagonal terms give the single-particle coupling between the two 
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levels. The real, symmetric matrix g aT describes the scattering of the atoms in 
each component off one another and the cross-scattering between components. The 
single-particle terms are trivial to treat, so the only new effect here comes from the 
cross scattering described by g l2 . 

The next step is to go to the interaction picture where the condensate mode is 
displaced to vacuum, just as in Eqs. (3.17). In this interaction picture, the field 
operators transform according to 

tMx)Z>(a,0(t)) = aUt)0«r(x,t) +4b(x) , (3.70) 

thus allowing an expansion in powers of 1/N 1 ^ 2 = l/|a|. We can write this transfor¬ 
mation more abstractly as 

£> T (a,0W) tp£>(a,0(t))^ = «I0W) + I > > (3.71) 

The interaction-picture Hamiltonian, as in Eq. (3.22), is given by 

'Hintit) = -ih ( ) + a (rb | 0(f)) - a* (0(f) | rj;) ) + D T (f) U(t) V(t) . 

(3.72) 

The time derivative of the condensate state is 

I 0W ) = “ (“1 I ^1 )) 0 I 1 ) + J t ("2 I 02 )) 0 I 2 )) • ( 3 -73) 

Putting all this together, we get the interaction-picture Hamiltonian to Bogoliubov 
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-i^+Ha+l \4 >t\ 2 )^Ao 


'Hint (t) = f ^2 a* 4>\ 

a 

+ 

T 

+ / (- ift s+^+E 

“I - ^ ^ fUjJ (TT Oi T (j). 

T 

+ J f + 5^^<7r|«r| 2 |0r| 2 ) 

“I - ^ ^ (ji^crT 9<JT^a4 ) a^r ( t , r^^ 


dx 

9ar\ a r\ \ *Pr I ) 
+ H.c.^ dx 


+ H.c. j j dx , 

CT,T / 


where the single-body translational Hamiltonians are 


H 2 2 

H = -V + V 

±J - cr o v 1 K a • 

2m„ 


(3.74a) 


(3.74b) 


(3.74c) 


(3.75) 


Just as for the case of a single component, we can neglect the c-number, mean- 
field-energy term. By requiring the term of order N 1 ^ 2 = |a| to vanish, we get a pair 
of coupled GP equations, 



ih Ft + h ’ + T, 


9<Tt\^7 



OL a 4>a + ^r«r 4>t = 0 • 

T 


(3.76) 


Notice that these are best thought of as coupled equations for the unnormalized wave 
functions, a 1 (j) 1 and a 2 4> 2 - It is often convenient to have the two GP equations written 
out separately as 


/ ^ \ 

ih— + idg^ J «i0i + hu 12 ct 2 (j) 2 = 0 , (3.77) 

\^Wt + ) “ 2< ^ + ^ 21 = 0 ’ (3-78) 
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where the GP Hamiltonians are 

-Hgp = Hi + hw u + 5 , n|tt 1 | J |</> 1 |~ + fji2 1Qy1 2 1021^ > (3.79) 

H^ = H 2 + hw 2 2 + ^72210^21 2 1021 2 + 92iWi\ 2 \^i\ 2 (3.80) 

(remember that oj 21 = UJ 12 and g 2 1 = g\ 2 )- It is also convenient to compactify the 
equations in terms of spinors relative to the two hyperhne levels so that we can take 
advantage of our bra-ket notation, 

where 

(hH hu 12 \ 

H gp = ^ = -^gp 11) (11 + -^gp I 2) (2 | + hu 12 11) (2 | + hu 2 i | 2) (11 . 

V^ 21 hwj 

(3.82) 

Recognizing that the spinor in Eq. (3.81) is the spinor representation of the state 
a | </>), we can write the coupled GP equations in the very compact form 

\ ih Qt + H zp) 1^) = °’ ( 3 - 83 ) 

where it is assumed, as our formalism requires, that a does not change in time. 


The Bogoliubov Hamiltonian governing the dynamics in the interaction picture is 
given by Eq. (3.74). In 4 x 4 matrix form, we have 

( |Mr)\ 


qj — _ 

rL bog r) 


((^l| (^2 I (*l>l| {A>1\) H bog 


A>2) 

iVl) 

VM4>/ 


where the symplectic-style matrix H hog takes the form 


(3.84) 


-ffgp + 9nK| 2 |di| 2 
U°21 + 921 a * a 2 0102 

Udl2 + 9l 2 aiOt 2 <^!02 

Hgp + 922 \ a 2 2 02 2 

9n a i 0i 

921 a l a 2 0102 

9l2 a l a 2 0102 

922 a 2 02 

9n ( a i) 2 (bi) 2 
g 12 a*ia 2 010 2 

92i a l a 2 0102 

922 ( a 2 ) 2 (^ 2 ) 

Hgp +9ul«i| i |di| 2 

b^>12 + 9l2 a l a 2 0102 

hui 2 i + 921 a l a 2 0102 

Hgp + 9221 a 2 1 1 02 1 
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(3.85) 


To get back to the compact spinor notation, we introduce, along with the matrix (3.82), 
two other matrices that operate in the spinor space defined by the hyperhne levels 
| 1) and | 2 ): 



(3.86) 


With these matrices, we have 



Notice that since $ is diagonal and G is real and symmetric, d'G'd 1 * and d>*Gd> are 


both Hermitian, and they are transposes and complex conjugates of one another; TGT 
and d>*Gd>* are both symmetric, and they are complex conjugates and Hermitian 
conjugates of one another. Using our total field operator and interpreting the 2x2 
submatrices as operators in the space of the internal levels, we can write the Bogoliubov 
Hamiltonian in the suggestive form, identical to that for a single component, 



(3.88) 


To eliminate phase diffusion in the condensate mode, we now introduce the 
auxiliary (nonHermitian) Hamiltonian T in exactly the same form it has in the 
single-component case [see Eq. (3.34)], 


r{t) = (A r-Nf + (c*4 (t) + - N)r ± + (a*a m - N)P ± , (3.89) 


where Af± = A f — and where the coefficient rj and the operator are defined 
in analogy to the single-component case, 


r)(t) = (0 |d>Gd>*| 0) = (0* |d>*Gd>| 0*) = (0* |<TG<r|0) = (0|d>Gd>|0*) (3.90a) 
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and 


J 7 1 



-a(xl> ± |$G$*| 0 ) 

7?a*a 0 + //era* - a*( <j> \ 1 1 |>) - a( ij) || 0) 



5^^crr|a<T| 2 |0<T| 2 QV0T^T± + H - C -j » 

rr t ' 


(3.91a) 

(3.91b) 

(3.91c) 


where | ) is the total held operator with the condensate mode excluded [see 

Eq. (3.67)]. Just as in the single-component case, rj is of order 1/N and J r ± is of 
order 1/N 1 ^ 2 . 


The transition to the interaction picture goes exactly as in the single-component 
case, yielding Eq. (3.39) at Bogoliubov order N°. Dropping the c-number term from 
the result, we have 

^int(^) = ( 2| a\ 2 a\a^ +a 2 a\a\ +(a* fa + aa\ + a J 7 } (3.92a) 

V f n \ |2 t , 2 t t , / *\2 \ 

— - I 2|a| + a a^a'^ + (a ) ) 

z y 7 (3.92b) 

+ ((M 2 a$, + (a*) 2 a 0 ) (0 | ^ ) + H.c.) 

+ (a*) 2 (rj; 1 | (f >*) ( (f>* |<TG<r | 0) (0 | ) + H.c.) 

7 (3.92c) 

- (|a| 2 (rH</>)(c/>|$G$*|rj>} 

+ (a*) 2 (r|d | (j )*) ( (f>* |<T 1 1 |>) + H.c.) . 


In symplectic form, we have 


^int (t) = g ' ((tl> | (4» T | ) F mtO) 



(3.93) 


Here the matrix of symplectic structure is given by 


F int(0 


1 \a\ 2 (Q$G<h*Q-<hG<h*) a 2 (Q<hG$Q* - ^ 

v (a*) 2 (g*<TG4>*g - <TG<r) |a| 2 (g*4>*G$g* - <TG<l>) y 


(3.94) 
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where Q(t) is the projector onto the single-particle subspace orthogonal to the 
condensate mode, as in Eq. (3.46). The modified (number-conserving) Bogoliubov 
Hamiltonian matrix assumes the form 


^ncb r, 


<4>| <^|) 


H 


neb 


i+) 


(3.95) 


with 


H = H 

''neb ''bog 


^"int — 


H gp + |a| 2 Q$G$*Q a Q&G&Q* 


a*) 2 Q*$*G$*Q H* + \a\ 2 Q*$*G$Q 


2 /~\* i.* / 


(3.96) 


The related demonstrations that at Bogoliubov order, if the condensate mode 
begins in a coherent state, it remains in a coherent state and that the auxiliary 
Hamiltonian tF(t) does not change the evolution in the iV-particle sector can be 
repeated word for word from the single-component case. There is, however, an 
important difference from the single-component case, which involves the orthogonal 
mode [ 4)) of Eq. (3.56). We get at this difference by first plugging into the number- 
conserving Bogoliubov Hamiltonian (3.95) the field-operator decomposition (3.65c); 
this separates out the condensate mode | (ft) and brings 7f ncb into the form (3.52a). By 
further separating out the mode | (ft) using the field-operator decomposition (3.65b), 
we arrive at 

n nch = aJa^l-Hgpl^) + (aj,a 0 (0|# gp |0) + 4( (ft |i7 gp | ) + H.c.) 

+ 4®4^ |(-ffg P + |a;| 2 $G$*)| (ft) + (a\,{(ft \(H gp + |aj 2< f ) G < f ) *)| 4>j_) + H.c.) 

+ ^(« 2 44(^ l $G$ l ft ) + « 2 4(4 $G$ I^jl ) + H.c.) 

+ ^ncbJL > 

(3.97) 

where 

^ncbJL = x I ((^jl | (tli |) H bog ( ] i (3.98) 

2 \\o) 

By the same calculation as for a single component, this Hamiltonian conserves a^. 
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We want to study the evolution of the mode | 0); for that purpose, we need only 
keep the parts of the Hamiltonian (3.97) that depend on a^ and aj,. Re-organizing 
the terms a bit, we have 


■Hncb = + H-C-) 

(3.99a) 

+ ( a\>a ( p( 0 \H gp \ ) + H.c.^ 

(3.99b) 

+ ( \ a 2 hJ,( 0 < hG < h* rpjL) + - a 2 a^(0 ij>j_ ) + H.c. j , 

(3.99c) 

where the term (3.99b) evolves the state 0) according to the GP equation, while 
the term (3.99c) expresses the coupling to modes orthogonal to both 0) and 0). 
The remaining term (3.99a) takes the form 

nj {oi 1 2 — f — , 2 id 2 —t —t , -2 10 / * \ 2 — - \ /Q inn \ 

H ss — ^ ^ ^ & (Ctf ) CLfpCLfpJ ^o.lUUciJ 

V ( id -f . —id *- \2 hl®l 
— 2 \ e aa 0 + e a a <f>) 2 ’ 

(3.100b) 

where 


? 7 (t) = (0|$G$*|0> = e~ 2ie (0 0*) = e 2ie (0* |$*G<r|0) 

= j J, 4 j (gn\<i>l\ "h ^22 1 02 1 2^ 12 |0 1 | |0 2 |)dx, 

(3.101a) 

(3.101b) 


with 9 = arg(aqa 2 /a 2 ). 


Presuming that aq and a 2 do not change in time—and, hence, that 9 is also 
constant in time—we can choose a e l6 = 0 to be real and positive, thus putting the 
Hamiltonian in the form (after discarding the irrelevant c-number that comes from 
operator ordering) 

?4s = (4 + ci^Y = fj/3 2 xl , (3.102) 

where x$ = (aj, + a^)/\/2 is the position quadrature corresponding to a^. The Hamil¬ 
tonian (3.102) produces shearing and squeezing in the direction of the momentum 
quadrature at a variable rate given by 2 fj{t)/3 2 . 

Another way to think about the Hamiltonian (3.100) is to recall that it arises, 
in the Bogoliubov approximation, from replacing a$ and a]j by a and a* in the 
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original Schrodinger-picture Hamiltonian. Restoring, in normal order, the creation 
and annihilation operators for the condensate mode to the Hamiltonian (3.100a) gives 
a Kerr-like interaction 

n ss = ^ ^244°0^0 + e 2 * e 44 a 0 a 0 + e 2l9 (3.103a) 
— ^ (e* e <4 a 0 + e ^ (4 a 0 + . (3.103b) 

The first term in Eq. (3.103a) comes from scattering of 0- and ^-particles off one 
another, the second term from scattering of two 0-particles into the 0-mode, and the 
last term from scattering of two 0-particles into the 0-mode. 


As an example of how the Hamiltonian (3.103) works, consider the situation where 
the condensate wavefunction 0 is an equal superposition of the two hyperhne levels, 
i.e., a ± = a 2 = l/-\/2, 6 = 0, and 


a <t> ( ai ° 2 ) ’ “0 ^2 


1 


a j 


= 71 ( 0l ~ ^ ’ (3.104) 

where Gq(a 2 ) is a shorthand for ai < p 1 (a 2t< j > )- By introducing the Schwinger operator 

0 'z — 2 (^iTi ^2^2) ^ (^0^0 ^0^0) 5 ( 3 . 105 ) 


we bring the Hamiltonian (3.103b) into the form 

Mss = 2fj(t)J 2 - | {a\a^ + 4 ^ 0 ) • (3.106) 

The J 2 term is the so-called one-axis-twisting [KU93] Hamiltonian and is widely used 
to generate spin squeezing in BECs [EGW + 08, GSH09, RBL + 10, GZN + 10]. The 
interplay of spatial and spin dynamics is considered in [SC00, SDCZ01, LTRS09]; 
Sinatra et al. [SWD + 11] showed that the amount of squeezing is bounded from above 
by the initial noncondensed fraction at finite temperature. 
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Chapter 4 

Bosonic Particle-Correlated States 


It would indeed be remarkable if nature fortified herself against further 
advances in knowledge behind the analytical difficulties of the many-body 
problem. 

- Max Born 


4.1 Introduction 

Quantum many-body problems are notoriously hard. This is partly because the 
Hilbert space becomes exponentially large with the number of particles N. As 
a consequence, one needs an exponentially large number of parameters merely to 
record an arbitrary state, not to say computing its time evolution. While exact 
solutions are often considered intractable, numerous approximate approaches have 
been proposed. A common trait of these approaches is to use an ansatz such that 
the number of parameters either does not depend on N or is proportional to N, 
e.g., the matrix-product state for spin chains [Whi92, Whi93, VPC04, DKSV04], the 
BCS wave function for superconductivity [Coo56, BCS57a, BCS57b], the Laughlin 
wavefunction for the fractional quantum Hall effects [Lau83] , and the Gross-Pitaevskii 
theory for BECs [Gro61, Pit61]. 
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The Gross-Pitaevskii theory, which uses the product ansatz, has precisely predicted 
many useful properties of Bose gases at ultra-low temperature. As particle-particle 
correlations become important, however, it begins to fail. To capture the quantum 
correlations, we propose a new set of states, which constitute a natural generalization 
of the product-state ansatz, 


p( x i , 


(0 x (0 WO . V (Q v (0 

i ^2 J * * * ? i J 1 5 J2 > 

r (0 „(0- 


7 ? s ,cr ( x i\ Yi ) <8> cr( x 2 , y 2 ) 


y W 

5 ,7 n 

(0 „(!)> 


) 


v ( 0 ')-p 

^' x n ) y n ) ' S 


tr (V 5 cr(x[°, yj°) ® a(x ( 2 l} , y ( 2 l> ) ® ■ • • ® cx(x^, y?) 


(4.1) 


where x^ = (x jA , x i>2 ,..., x jV ) and = (y i4 , y ii2 , • • • • Yjj) denotes the coordi¬ 
nates of l particles, cr is an arbitrary state (density matrix) of the l particles, 1 and Vg 
is the projection operator onto the symmetric subspace of all the N = l x n particles, 


'Ps | ^1, "02, ■ ■ ■, i’N ) = ^ | Vhr(l)> ^(2), ■ • • > VV(iv) > , (4-2) 

ttSSjv 

where the sum is over the permutations n in the symmetric group S N . The state (4.1) 
is derived by symmetrizing the n-fold tensor product of the /-particle state a; we call 
the resulting state a Bosonic Particle-Correlated State (BPCS). 2 As a consequence of 
symmetrization, the quantum correlations existing in the /-boson state cr “spread out” 
to any subset of the / x n bosons. The parameter space of the BPCS does not grow 
with n; it equals to that of the bosonic states for / particles. 


In this dissertation, I pay most attention to the special case that cr is pure and 
/ = 2. When / = 2, PCS can also be read as Pair-Correlated State. 


In this chapter, I show that BPCS is a many-body ansatz which is both quan- 
tumly correlated and computationally efficient. One advantage of BPCS is that 
it can represent quantum states with or without Off-Diagonal Long Range Order 
(ODLRO) [Yan62], For example, both the superconducting and the Mott insulating 

1 The /-particle states cr can be restricted to symmetrized states, or they can be left 

arbitrary, with Vg taking care of the symmetrization when the BPCS state is constructed. 

2 

For the fermion case, we can simply substitute the anti-symmetrizing operator V^ for the 
symmetrizing operator Vg. The resulting state can be called a Fermionic Particle-Correlated 
State (FPCS). 
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phases in the Bose-Hubbard model can be described by the BPCS ansatz. An inter¬ 
esting question is whether the BPCS ansatz can faithfully interpolate the two phases. 
In chapter 4, we will show that the answer is yes even if we restrict ourselves to the 
case 1 — 2] there we compare both the ground state and the dynamics predicted by 
BPCS to fully numerical results. 

4.2 Historical Remarks 

A central topic in condensed matter physics is to study the particle-particle correlations 
existing in many-body wave functions. Here I briefly review some existing approaches 
that are used to capture the particle correlations in ultracold bosonic gases, with an 
emphasis on the relations to the BPCS ansatz. 

One of the first and most influential approaches to BECs that take interparticle 
correlations into consideration is the Bogoliubov approximation [Bog47, Fet72, Hua87]. 
Although the Bogoliubov approach has made many precise predictions, it is a pertur¬ 
bative method; i.e., the depletion of the condensate must be small for it to work. I 
will argue that the nonperturbative BPCS recovers the number-conserving Bogoliubov 
approximation in the limit of small depletion. 

Inspired by the BCS wavefunction proposed by Bardeen, Cooper, and Schrieffer 
for superconductivity [BCS57a, BCS57b], Valatin and Butler [VB58] introduced a 
similar pairing wavefunction for bosons, 



(4.3) 


where N is a normalization factor. This state, with a quadratic form of the creation 
operators in the exponential, is very different from a coherent state. The coherent state 


has a Poissonian number distribution which is peaked around N = |a| 2 , whereas the 


Valatin-Butler state (4.3) satisfies an exponential number distribution. Consequently, 
it might not be suitable to describe situations where the total number of particles is 
conserved. One easy way to see the exponential number distribution is by setting 
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A j = 0 for all j ^ 0, which gives 


exp (A 0 aJaJ) | vac) 

(4.4a) 

-j 00 \Tl 

1 A 0 / f\2n , , 

> —r {%) vac) 

v n=0 

(4.4b) 

-j 00 \ n 

- 7 = V -T V(2n)! |2 n)„® vac) ± 

(4.4c) 

1 ^ (2A 0 ) n . v . . 

.— 2n) 0 <g) vac) ± , 

v n=0 

(4.4d) 


where we use Stirling’s formula in the last step. It turns out this exponential 
distribution is more general; it is valid even when more than one A j is nonzero. 3 
Unlike the pairing wavefunction for fermions, which is always normalizable, the 
Valatin-Butler wavefunction (4.3) can only be normalized when |Aj| < 1/2 for all 
j G {0,1,...}. The Valatin-Butler wavefunction has been used to investigate the 
transition from a single condensate to a multicondensate [CM67, EI69, NSJ82], 


A number-conserving version of the Valatin-Butler wavefunction was introduced 
by Leggett [LegOl, Leg03], 

1 / t t . ^ . + + \ 


^legg ) 


a/M V 


( a Wo + 2 AfcCl fc a 


t 

—k 


vac) , 


(4.5) 


where the normalization factor is not exact. This state is a special case of the BPCS 
with a pure and 1 = 2, but a systematic treatment seems to be lacking in the literature. 
I will discuss the single-particle and two-particle reduced density matrices of this 
IV-particle state in the following sections; a generalized Gross-Pitaevskii equation 
is derived using this ansatz in the next chapter. It is worth mentioning that the 
ground state of a spin-1 Bose gas with an antiferromagnetic interaction takes a 
form similar to Eq. (4.5) [HY00]. The fermion version of this state was also used 
by Leggett et al. to treat the BEC-BCS crossover problem [Leg80, LZ12] and by 
others [CLT03, Law05, CBMD08, TBM13] to discuss the composite boson problem. 


3 In the large n limit, the normalization factor of the 2n-particle sector scales as (2"n! A”) 2 
[see Eq. (4.46)]. By replacing \J (2n)! in Eq. (4.4c) with 2 ra n!, we have the same exponential 
distribution as in Eq. (4.4d). 
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It is difficult to do any calculation with the form (4.1), because of the need to do an 
explicit symmetrization. By going to the second-quantized picture, symmetrization is 
done automatically. For a pure PCS with l = 2, the PCS is specified by a two-boson 
wavefunction !7 r(2 ' ) (x 1 , x 2 ), and the PCS is given by 

^ P cs(xi,x 2 ,...,x 2n ) oc x 2 )^ ( 2 ) (x 3 , x 4 ) ■■■& (2) (x 2n _ 1 , x 2n )) . (4.6) 

Such a PCS can be regarded as constructed by a mapping of the two-boson Hilbert 
space into a submanifold of the Hilbert space of 2 n bosons. Note that the two-boson 
wavefunction always has a Schmidt decomposition of the form [PY01], 

V 

^ (2) (x 1; x 2 ) = J^Aj-^XiN^Xa) ’ ( 4 - 7 ) 

3=1 

where v is the Schmidt rank and the A^s are the square roots of the Schmidt coeffi¬ 
cients (Xq=i A 2 = 1). The coincidence of the Schmidt bases of the two bosons is a 
consequence of the symmetry of the wavefunction. 


It is instructive to perform the Schmidt decomposition in the second-quantized 
picture, where an arbitrary two-boson state takes the form 

I ^ (2) ) = ~Jd A P h k h ) I vac ) • ( 4 - 8 ) 

V Z j,k =1 

Here fet is the creation operator of the jth single particle state, and A jk = A k j is a 
symmetric 4 matrix, normalized according to tr(AA^) = 1 to make | ) normalized. 

The Autonne-Takagi factorization theorem [see Corollary 4.4.4 (c) of [HJ13], or App. E 
of this dissertation] says that any complex symmetric matrix A can be diagonalized 
by a unitary matrix U in the following way, 

UAU T = diag (A x , A 2 ,..., A„) , (4.9) 


where the A^s are real and positive and normalized as Xq=i = 1- The XjS are the 
singular values of A, and the diagonalization (4.9) is a special case of the singular-value 

4 We can always make A symmetric by redefining A —> (A + A T )/2, without changing 


t^ (2) ). 



Chapter 4. Bosonic Particle-Correlated States 


decomposition, specialized to symmetric matrices. We adopt the following convention 
throughout this dissertation, 


Ai > A 2 > ••• > A, . (4.10) 

By introducing a new set of creation operators = Ylk=i Ujk we l iave 

l^ (2) > = “4 EE®])" | vac) ; (4.11) 

j = i 

the corresponding wave function ( x l5 x 2 | ) = l^ r ^ 2 ' ) (x 1 , x 2 ) has the Schmidt 

form (4.7), with ipj (x) being the single-particle wavefunction that goes with the 
annihilation operator a 3 , i.e., 

dj — f i/j*j (x) rjj(x) dx . (4.12) 


We now define the pair creation operator 

•^XX“1) 2 = 

3 = 1 

Using 


J i^ (2) (} 


■1 5 


2 )rl) T (xi)^ t (x 2 )dx 1 dx 2 


(4.13) 


4' t (x 1 )---'h t (x w ) | vac ) = VhP.V s \ Xij.-.jXjv) (4.14a) 

— l X 7r(l)’ • • • > X 7T(iV) ) > (4.14b) 

V ■ ttGS n 

where the sum is over the permutations 7r in the symmetric group S N , we have 
(^4 t ) n | vac ) 

= vOV! J | x„..., x 2 „ > P s (^ <2 > (x,, x 2 )... i' 12 ’ (x 2 „_„ x 2 „)) dx,... <lx 2 „ , 

(4.15) 


where 


7> s (V 2) ( x i, x 2 ) ^ (2) (x 2n _!, x 2n ) 


IV! 


E ^ (2) ( X 7r(l)) X 7T , 2 ,)--^ ,2, (x i(2 „_ 1| ,x i(2 „ i ). 


TTSSa 


(4.16) 
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Equation (4.15) is the relation between the first- and second-quantized pictures; it 
can be written in the equivalent form 

( x i) ■ ■ ■;x 2n |(«4 t ) n | vac ) =P 5 (V 2) ( Xl , x 2 ) ■■•lP r(2) (x 2n _ 1 , x 2n )) , (4.17) 

which states that one gets the PCS state (4.6) by applying the pair creation operator 
n times to the vacuum. 


The state in Eq. (4.17) is not normalized. The properly normalized PCS state is 
given by 


S'pc) = T (V)" I vac) , 


(4.18) 


where N is a normalization factor that plays an important role in our consideration 
of PCS states and n = N/2 is the number of pairs. We abandon the normalization 
restriction on the Schmidt coefficients from now on, requiring only that the XjS be 
real and positive, since an overall scaling of the XjS is automatically absorbed into N. 


The form (4.18) is convenient for calculations, but let us first build some intuition 
by considering the state decompositions in the first-quantized picture. The relative- 
state decomposition of the particle x 4 relative to all the other particles is 

^pcs OC PsfV^Xr, X 2 ) ^ (2 ) (x 3 , X 4 ) X 2n)) (4.19a) 

« X J Vh'( X l) ^s(^'( X 2) ^ (2 ) ( x 3, x r) ■ ■ ■ ^ (2) ( x 2n-l, x 2n)) ■ (4.19b) 

3 

This is a Schmidt decomposition, and the Schmidt basis of the particle x 4 consists of 
all the single-particle wavefunctions ^ J (x 1 ). The Schmidt coefficients, however, are 
not given by the A^s, because the norms of the relative states of the other particles 
are different for different values of j. 


More interestingly, we have the following relative-state decomposition for particles 
x, and x 2 , 

^ P cs oc fl^Xi, x 2 )P 5 ^ (2) (x 3 , x 4 )---^ (2) (x 2n _ 1 , x 2n )j (4.20a) 

+ ( N ~ 2 ) X i Xk ^'( x i)^fc( x 2 ) 

jX 

X 'Psi'I’jM'l’kM ^ ( 2 ) ( x 5 , x e) • ■ ■ ^ < 2 ) ( x 2 n-l, x 2 n)) ■ 


(4.20b) 
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What this shows it that the reduced (marginal) two-particle state consists of two 
terms: In the first term, the two particles x 1 and x 2 can be perfectly correlated, 
whereas in the second term, they are only partially correlated. To determine the 
pairwise quantum correlations in the PCS, we need to find the two-particle reduced 
density matrices (2RDMs), and to do that, we will find it useful to investigate the 
normalization factor N. Before getting to that, however, we detour into showing how 
the Bogoliubov approximation arises from the PCS. 


4.4 Relation to Bogoliubov’s Wavefunction 

When the Schmidt rank of the two-particle state (4.11) is one, the PCS (4.18) 
is a product state (this is the Gross-Pitaevskii ansatz). Here, I show that the 
PCS also reproduces the particle-number-conserving (^conserving) Bogoliubov 
approximation [GA59, Gar97, CD98] to the BEC ground state when there is one 
dominant and other minor Schmidt coefficients. This result validates using PCS for 
BECs for small depletion, and it also provides a different way to write the state of an 
N -conserving Bogoliubov approximation. 

Consider a Bose-Einstein condensate of N particles in a trapping potential. The 
Gross-Pitaevskii (GP) ground state takes the form 

l^gpgs) = I N)o® |vac)j_ , (4.21) 

where the subscript _L denotes all the modes that are orthogonal to the condensate 
mode. The essence of the IV-conserving Bogoliubov approximation is to perturb 
about the state (4.21) by introducing the operators 

a] = ata 0 / VN for j — 1, 2,..., v — 1, (4.22) 

where a 0 is the annihilation operator of the condensate mode and a]- is the creation 
operator of the jth orthogonal mode. Note that even though the operators oj and 
hj are not exactly creation and annihilation operators, they satisfy the canonical 
commutation relations approximately in the limit of large N and small depletion 


Chapter 4. Bosonic Particle-Correlated States 


61 


(small number of particles not in the condensate mode). For a condensate with small 
depletion, a 0 is of order y/N and hjS are of order 1; the modification to the total 
energy in the Bogoliubov approximation is of order 1. 

The IV-conserving Bogoliubov Hamiltonian, with the mean field removed, is a 
quadratic function of a] and cq, 

u-i i 

Tdnch — Mj k Oja k + - (^Mj k a^a\ + H.c.^ , (4.23) 

j,k=l 

where = M and (Ad') 1 = M 1 . This Hamiltonian can always be diagonalized by a 
Bogoliubov transformation 

V -1 

B 1 H nch B = ^e k a\a k , (4.24) 

fc=i 

where B is a Gaussian unitary inducing a symplectic transformation on the operators 
a]- and aj for j — 1,2,, u— 1. The Bogoliubov ground state of the Hamiltonian (4.23) 
thus takes the form 


^b g s) = B | N ) 0 ® | vac ) ± 


u-l 


U ( n s kilk, a k ) J v f I N ) 0 <g) I vac ) ± 
u ( n S k {lk, a k )\ W\N) 0 ®\ vac) L , 


k =1 


(4.25a) 

(4.25b) 

(4.25c) 


where we use the Bloch-Messiah reduction theorem [BM62] to decompose the Gaussian 
unitary B into a multiport beam splitter V , followed by a set of single-mode squeezers 

1/ -2 ~f 2. 

5 fc ( 7 fc , a k ) = e 2 l7fe “ fc lkCLk ’ with 7 k real, followed by yet another multiport beam splitter 
U. Also note that W and IT have no effect on the GP ground state | N) 0 <g) | vac)j_; 
thus we can use either of them (or neither). 


The action of the multiport beam splitter U is given by 

V- 1 

U a.j IP = a k U k j = bj = a^bj /x/iV , 

k =1 


(4.26) 
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where U is some unitary matrix that specifies the multiport splitter. Note that bj is 
the annihilation operator for a new set of modes orthogonal to the condensate mode. 
The state (4.25) thus takes the form 


V -1 


' k =1 ' 

v—1 \ 

n s k(lk,h)j I N) 0 ® I 

k =1 

l /—1 

n 


vac 


vac 


1 / tanh 7 fc ~ t 2 


fc= i \/cosh 7 fc 

1 

rifc=i c ° sh 7fc 


exp 


exp 


2 

i/-i 

£ 

k =1 


tanh 7 fc If2 



(4.27a) 


(4.27b) 

0 ® I vac ) ± 

(4.27c) 

N) 0 <g> | vac ) ± . 

(4.27d) 


Here we use the “quasi-normal-ordered” factored form of the squeeze operator [Per77, 
Hol79]: 

S(7, a) = exp ta1 ^ 7 a f 2 ) (cosh 7) “ exp a 2 ) . (4.28) 

Note that we can always make 7*. negative, by redefining the phase of b k , so that the 
coefficient | tanh j k in Eq. (4.27d) is positive. 


On the other hand, when there is one dominant Schmidt coefficient A 0 , the 
PCS (4.18) of N = 2n particles takes the form 



(4.29a) 

(4.29b) 

(4.29c) 


where the approximations are good when N is large and the depletion from the 
condensate mode is small. Comparing Eqs. (4.27d) and (4.29c), we find that they 
can be made the same by choosing 


V -1 

N = N\ A ( f cosh 7 fc , 

k =1 


(4.30) 


A j /A 0 = — tanh 7^ , for j — 1, 2 ,..., v — 1 . 


(4.31) 
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Hence, as promised, the PCS (4.18) encompasses the Bogoliubov approximation to 
the BEC ground state. 


Another way to prove the same result is by noticing that 




pcs 


) ~ v * exp (w £ A 


v-l 


k =1 


a 


vac 


(4.32) 


where V N is the projection operator onto the IV-particle sector and | a ) 0 is a coherent 
state for the condensate mode with a = \fN. This is nothing but the extended 
catalytic state (3.5) that we have introduced for the number-conserving Bogoliubov 
approximation. 


4.5 The Normalization Factor 


The importance of the normalization factor N to PCS states is analogous to the 
utility of the partition function in statistical physics. By taking derivatives of the 
normalization factor with respect to the Schmidt coefficients Xj, one can calculate 
the reduced density matrices (RDMs) in the Schmidt basis, and these, in turn, give 
all the physical observables. 

In the second-quantized picture, the normalization factor N introduced in Eq. (4.18) 
takes the form 


N = (vac | A n («4 f ) n | vac) 


77 J ( vac \A n \a )( a | (^)"| vac ) d 2 aq • • • d 2 a v 



(4.33a) 

(4.33b) 

(4.33c) 


where we insert a complete basis of coherent states. Expanding the monomial in 
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Eq. (4.33c), we have 




n\ 
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rah 


Ui\u 2 \ ■■■u v \ 


n 

3 = 1 


(2uj)\ A 


2Ui 


Uj \ 


(4.34a) 


(4.34b) 


where Uj (vj) are non-negative integers satisfying Y^j=i u j = Y^'j=i v j = n - Although 
the sum in Eq. (4.34b) appears to requires exponential time to evaluate, it can be 
evaluated in polynomial time by using an iterative algorithm. It is still, however, 
computationally demanding for large N, in addition to not being intuitive. 


To make progress in interpreting and evaluating N, we take what might be 
construed as a backward step by writing N in a different integral form. To do so, we 
use (2 Uj)\/uj\ = 2 Uj (2uj — 1)!! to write 

<-> 

and then use this to put Eq. (4.34b) in the form 

r\TL l /»00 o / U \ U 

N= pmf L e ~ w ' 72 (g A? *0 dVv '■ ^ - (436) 

We use the above expression to evaluate the normalization factor in the large-A r 
limit in the next section. For now, however, we employ it to derive several exact 
expressions. 



The first use of Eq. (4.36) is to derive a generating function, 


2 2 n i on 

n! o 


N = 


(v^T dr 

d 


-\y\ 72, 


exp (b r ) d yi--- d v. 


3 = 1 


T —0 


= 2 2n n\ 


dr r 


Ov/l-rA 


T—0 


(4.37a) 

(4.37b) 


An equivalent way to obtain the generating function (4.37b) is to evaluate the 
quantity (vac \e^ F ' A e y d FA |vac), which can be calculated using the “quasi-normal- 
ordered” factored form of the squeeze operator found in Eq. (4.28). 
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The second use is to derive expressions for the diagonal elements of RDMs. 5 As 
a first example, we are able to represent, using Wick’s theorem, the jth diagonal 
element of the 1RDM p ^ with the normalization factor, 



— ( vac | A n OjCij ( > A t ) n | vac ) 

—^ vac | A! 1 - 1 a 2 j (.4.^)"| vac ) + ( vac \A n (a]) 2 (M^)" X | vac ) j 

N ( ( vac I §X7 DTI > + < vac S tiL| vac )) 

\j ON 
'NdX'j ' 


(4.38a) 

(4.38b) 

(4.38c) 

(4.38d) 


The two terms in Eq. (4.38b), which correspond to contracting aj and aj with the pair 
annihilation and creation operators, are equal. In addition, by using Wick’s theorem, 
it is not hard to prove that all the off-diagonal elements of p l> in the Schmidt basis 
are zero and, therefore, the normalization factor and its first derivative determine the 
1RDM. 


More generally, we have the following result for the diagonal elements of the 
g-particle RDM p^\ 


ph-pjy-p = d ( vac I* 4 " a l • • • a l% ■ ■ ■ a h ( 4 TI vac > 


N 

A?! ' ' ' ^jq 


d q N 


N 


dXj ■ ■■d\~ 

J 1 Jq 


(4.39a) 

(4.39b) 


This result can be proved by mathematical induction. We already have that Eq. (4.39b) 
holds for q = 1, so to show that it holds for all positive integers q is that if it holds 
for q, it is satisfied for q + 1. The inductive hypothesis is thus that 


<9 9 N 


DX h • • • dXj 


vac 


\A n 


\?i 


Jq 


Jl 


CL » CL si 
Jq Jq 


*Jl 


DTI vac ) . (4.40) 


’The (/-particle RDM is normalized to N\/(N — q)l = N(N —!)••• (A — q + 1) in 
this dissertation unless stated otherwise; normalized in this way, the RDMs are equal to 
correlation functions. 
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By taking derivatives with respect to of both sides of Eq. (4.40), we have 
<+ +1 N 1 


<9 A a • • • dX q ; A,- • • • Xa 

3 1 3 q -\-1 3 1 3 q 


/ I An— 1 2 t t 
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3 1 3q 
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x | ( vac \A n 1 a jq+i a) i ---a)a jq ---a h 


a i q+ 1 (^ t ) n |vac)+c.c.) 


\'l ' ' ' A?' 9 +1 


vac (A" at ■ ■ ■ at- 


-7l Jij+l -?9+l 




(4.41b) 

(4.41c) 


which is the required result. We note that for q > 1, the g-particle RDM is generally 
not diagonalized in the Schmidt basis. In Sec. 4.9, we show how to construct the 
entire g-particle RDM using only the diagonal elements calculated by the above 
method. 


Although we can avail ourselves of the power of Wick’s theorem to derive the results 
of this section, we can demonstrate the same results using only the commutators 
[cij, (A r ) n ] = 2n\jAj(A Jf ) n ~ 1 and [M™,aj-] = 2nXjA n ~ 1 aj, which imply that 

dj(A^) n \ vac) = 2n\jAj(A^) n ~ 1 | vac) , (4.42a) 

(vac | A n Oj = 2n\j(v&c \A n ~ 1 aj . (4.42b) 


4.6 The Large-iV Limit 

Often there are thousands to millions of atoms in a BEC, and it is sufficient to 
work with results that are valid in the large -N limit. In this section, we discuss 
how to derive an asymptotic form of the normalization factor for large N. To 
get the desired analytical results, terms that are of order 1/N smaller than the 
leading terms are neglected. This means that the following results do not include the 
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Bogoliubov approximation ; 6 instead we should think of them as a generalization of 
the Gross-Pitaevskii approximation. 
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Figure 4.1: The three eigenvalues of the single-particle RDM (normalized to 
unity) plotted as a function of the number of particles N for the case Ai = 0.5, 
A 2 = 0.3, A 3 = 0 . 2 . As N gets large, the eigenvalue corresponding to the biggest 
A approaches one while the other eigenvalues become negligible. 


Before the analytical calculation, let us look first at some numerical examples to 
get some intuition. In Fig. 4.1, the eigenvalues of the 1RDM (normalized to unity 
instead of to N) are plotted as a function of the number of particles N for A^ = 0.5, 
A 2 = 0.3, and A 3 = 0.2. For N = 2, the eigenvalues of the 1RDM are, of course, equal 
to the XjS, but as N gets larger, the eigenvalues become further apart. Eventually, 
the biggest eigenvalue approaches one, leaving the other two eigenvalues negligible; 
thus, as far as the 1RDM can tell, the PCS becomes an uninteresting product state 
for large N. 

In the second numerical example, plotted in Fig. 4.2, we consider the situation 
where there are two XjS that are nearly degenerate, Ai = 0.41 and A 2 = 0.39, and a 

6 The depletion predicted by the Bogoliubov approximation is of order 1, i.e., it modifies 
the 1RDM by order 1, which is 1/N times smaller than the 1RDM itself. 
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third smaller value, A 3 = 0 . 2 . When N is of the order 1 /(Ax — A 2 ) = 40 or larger, the 
third eigenvalue has died out, and only the two biggest eigenvalues play much of a 
role in determining the 1RDM. Our numerics suggest that in the large -N limit, only 
those AjS that are within order 1/N of A x survive. 



Figure 4.2: The three eigenvalues of the single-particle RDM (normalized to 
unity) plotted as a function of the number of particles N for the case A^ = 0.41, 
A 2 = 0.39, A 3 = 0 . 2 . The solid lines are the eigenvalues calculated by setting 
A 3 = 0 and keeping X 1 and A 2 unchanged except for renormalizing; they conform 
pretty well with the other results for large N. 


Generally, we speculate that one only need to keep those A^s that are within 1/N 
of A 1; the largest eigenvalue given our ordering convention (4.10); the other A^s can 
be omitted without affecting the PCS—i.e., the relevant low-order RDMs are not 
affected. This speculation is already supported by the above numerical results, and I 
will argue further for it based on the analytical results for the large N limit. The 
important A^s, being very close to each other, can be rescaled (i.e., they are no longer 
normalized to one) and parameterized as 

A^ = 1 + — or A,- ~ 1 + -4- , (4.43) 

J n An 

where the SjS are real parameters of order unity. Because of the rescaling, all the A_,s 
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are very close to 1, and their differences are of order 1/N. Putting Eq. (4.43) into 
Eq. (4.36), we manipulate the normalization factor N through the following sequence 


of steps: 


2 n n\ 


= __ . 

{V2^Y J- 


2 n n\ 


OO 

oo 


0 -\y\ / 2 


D -\y\ / 2 


Y^V 2 j ) d '!h ''' Yh 

3=1 


J-o 

2 n n\ r 

yMy ]- c 

2 n n\ r 


Y l 1 + Ct) y 3 ) dyi ''' d '!C 


3 = 1 


0 -W\ / 2 uvi 2n 


,-\y\ / 2 \rt\ 2n , 


1 + 


n \y\ 


Y s j y 2 j) dyi '" dy v 


(V^y J- 

2 n n\ f 

(V^k y Jo 


oo 

oo 


exp ( —j Y s 3 y2 3 ) d Vi'-- dy . 


-r /2 2n+v-l 7 

e ' r dr 


4 n n\ f v 

^ r (" + 2 


l»l 

r v 

/ ex p ( Y s 3 y f) dn 

Ay\=i v j= 1 J 

r v 

/ ex p (Y s .) y j) dn > 

J\y\=l ' 


(4.44a) 

(4.44b) 

(4.44c) 

(4.44d) 

(4.44e) 

(4.44f) 


where dCl denotes the area element on the unit [u — l)-dimensional sphere \y\ — 1. 
The only approximation here is to replace, in Eq. (4.44d), the power function by the 
exponential function. For each low-degree monomial of SjS the error in its expansion 
coefficient as a result of this replacement is of order l/n; such error only becomes 
substantial when the degree of the monomial approaches n. This is an excellent- 
approximation for our purpose of calculating low-order RDMs, because the high-degree 
monomials only affect high-order RDMs. 


Denoting the Gaussian integral in Eq. (4.44f) by 
1 


T(«) - 


2 tR /2 

1 


27T 


v/2 


/ ex P[Y s i y j) dfl 

y \y\=i V j=1 J 

/ OO f V \ 

ti{\y\ - l) exp ^ Y s 3 y jJ dy i ' 


•• dy u , 


(4.45a) 

(4.45b) 


the normalization factor takes the form 


N i> ~4”ri!r(r l +0T(s) . 


(4.46) 
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The greatest significance of this expression is that the dependences on n and on s 
(or A) factorize. According to Eq. (4.39b), the diagonal elements of the RDMs can 
now be expressed approximately, with errors of order 1/N, as 


V-A ^T(A) 


T(A) 9X h ■ ■ 

■ dX j 

Jq 

(InYXl ...Xl 

d q r(x) 

r(A) ds h ■ ■ ■ ds jq 

(2 n) q d q T(s) 



T(s) ds h ---ds jq ‘ 


(4.47a) 

(4.47b) 

(4.47c) 


The function Y(s) determines the PCS in the large -N limit, with all the ?r-dependence 
removed. Relative to exact expressions like Eq. (4.34b), the complexity of evaluating 
T(s) is dramatically reduced because of the removal of the n-dependence. 


As an example, consider the fcth eigenvalue of the 1RDM, which in the large- N 
limit takes the form 

Pkk = f y 2 k eX P ( it, S i y t) dn ■ ( 4 - 48 ) 

2n 7 Y(s) J\y\=i \ j=1 / 

Consider the case where X u is substantially less than 1 and all the other A^-s are close 
to 1 . In this situation s u is a negative number with large magnitude. Because of the 
suppression of the exponential function in Eq. (4.48), the magnitude of y u must be 
very small to contribute to the integral, which tells us that fi/l <C 2 n. On the other 
hand, we have Xq=i Vj — 1 f° r ^ ie other dimensions, so for j = 1,2,..., v — 1 can 
be calculated by neglecting the last dimension y u . In effect, the integral in Eq. (4.48) 
is reduced to an integral over a hypersphere of one less dimension. This argument 
can be easily generalized to higher-order RDMs, and it confirms our speculation that 
we need only keep those A^s that are within 1/N of X 1 . 

Although we have already made life easier by introducing T (s) , it is still a difficult 
task to evaluate the Gaussian integral (4.45a) over the hypersphere. Fortunately, we 
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can reduce the expression for T (s) to a single-variable integral. To do so, notice that 

/»oo 9 poo 

/ x‘ ,/2 - 1 e- rx r( x i)d x =- 

JO 


2 V ! 2 


v—l —TV /2 op I ' w \ 7 

r e I — dr 


(4.49a) 






r — s,- 


(4.49b) 

(4.49c) 


where we do the substitution x — r 2 /2 in Eq. (4.49a) and where t > s 1 for convergence 
(s x is the largest of the Sj s). Because Eq. (4.49c) is the Laplace transformation of the 
function x^ 2 ~ X T(ys), we have 

1 


T ( X s)=x 1 -' I2 C- 1 [ ] f 


4/ 


T — S 


l—u/2 /»A+ioo / u -i 


dr , 


(4.50a) 

(4.50b) 


where T 1 stands for the inverse Laplace transformation and the real parameter 
A > s-l for convergence. We have thus succeeded in reducing the high-dimensional 
integral (4.45a) to the one-dimensional integral (4.50b). 

For numerical calculations, one might find a straightforward series expansion of 
the function Y(s) to be useful: 


n 

Z7T J\y\=lj = x 


m dn 


E 


m i m 2 m,. 

P 1 O ^ O 1 / 

*1 *2 b u 


2ir v I 2 ^ 


(4.51a) 

y^y^-'-yl^dn. (4.51b) 




The integral in Eq. (4.51b) can be manipulated by a change of variables into the form 


2 mi 

Vi 


=i 


/ LXJ 

<S(lz7l - ■ • • yl mv dy 1 ■■■ dy v 

-oo 

= 2 J 5 ^ Zj — 1 ^ z™ 1 x l 2 ---z™ v ~ x l 2 dz\- 

= 2B[m x + l/2,...,m I/ + 1/2) , 


(4.52a) 
■ dz v (4.52b) 


(4.52c) 
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where 


, N UU r K- + V 2 ) 

B(m 1 + l/2,...,m„ + l/2) = iii f i__l__ (4.53a) 

= ~L m (4.53b) 

2 I (m + v/2) 

is the multivariable Beta function, with m = [notice that (—1)!! = 1]. 

Putting this back into Eq. (4.51b) gives 


T(a) = - 

7T 


v/2 


oo 




3 1 *2 




mil m 2 \ ■ ■ • mj 


B(m 1 + 1/2,, ?n„ + 1/2) . (4.54) 


For interested readers we also show here how to represent the function T (ys ) as 
a convolution. Note that 




r - Sj V7T 


X - 1/2 dx 


1 f°° 

=- 7 = / e rx G s .(x)dx , 


7T 


where 

for x > 0 , 

0 , for x < 0 . 

Putting Eq. (4.55b) into Eq. (4.49c), we have 


G,,(x) = 


V /2_1 e" r * 


T(x?)* = n(^/ 


v 1 2 / 
7T 7 JO 


/ e TX G s .{x)dx 

do 2 

e- T *(G» 1 *G Sl *-..*Gj(x)<(x, 


where * stands for the convolution, 


(G..»Gj(x)= / G,.(x-x')G, t (x')dx' ■ 


Doing the inverse Laplace transformation, we have 


%“) = 


X 


l-v/2 


IT 


v/2 


(4.55a) 

(4.55b) 

(4.56) 

(4.57a) 

(4.57b) 

(4.58) 


(G ai *G. a *---*Gj( X ). 


(4.59) 
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We now have four representations of Y(s), expressed in Eqs. (4.45a), (4.50b), 
(4.54), and (4.59). All of these turn out to be useful, and we use whichever is most 
convenient. 


We turn now to an exploration of relations among the RDMs that can be derived 
and expressed through the function Y(s). First, from the definition (4.45a) or from 
the Laplace transform (4.50b), we have 

T(s + 51) = e s T(s) , (4.60) 

where 1 = (1,1,..., 1) T and 5 is a c-number. The only effect of adding a constant 5 
to all the SjS, which changes the normalization of the A^s by v5/n, is to change Y(s) 
and, hence, N ?n by multiplying by a factor e 6 . This trivial fact implies that 
ar(s) _ 0T(s + Sl) 


£ 

k =1 


dsi 


05 


6=0 


= A?). 


(4.61) 


which applied to Eq. (4.47c), confirms the normalization condition for the 1RDM: 


Eh 1 ,? 


2 n J^OT(s) 

Y.-M = 2n 


T(s) ^ 0s K 

k =l V° ) k =l 

Equation (4.61) can be generalized to 


(4.62) 


£ 


0 d q T (s) _ d q T (s) 


^ ds k ds~ ■ ■ ■ ds~ Os a ■ ■ ■ ds~ 

k=l K Jq -7i Jq 

which corresponds to the following condition for the higher RDMs: 


V o {q+1) 

rjl-jqkJi—, 


(2 n 


, 9+1 


£ 


a tfUs) 


k =1 


jqk T («) ^ ds k ds h ■ ■ ■ ds jq 

(2n) q+1 <9 9 Y(s) 

Y(s) ds 3] ■ ■ • ds 3q 

— 2n Pj v ..j jyj ■ 


(4.63) 

(4.64a) 

(4.64b) 

(4.64c) 


Notice that we have the factor 2 n, instead of 2 n — q, in Eq (4.64c); this is because of 
the approximations we have used, which should be fine for large n. 

More powerful relations of the RDMs can be derived using the Laplace form (4.50b), 


Y (s) = — f 
K J 2m A 


A+zoo 


A—zoo 


n 




dr . 


(4.65) 
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For sj 7 ^ s k , we find 


<9 2 T(s) 


^A+ioo 


dsj ds k 2ni J A _ ioo 4 (r - sj)(t - s k ) \ ^ y/r - s. 


n 


dr 


r*A+«00 


2™ J A -ioo 4 (Sj - S k ) \T- Sj T S k ) V^ji V T ~ S r 


n 


2 (Sj - S k ) V ds 


dr (s) or(s) 


ds. 


(4.66a) 
dr (4.66b) 
(4.66c) 


Equations (4.47c) and (4.66c) together give the following relation between the single 
and two-particle RDMs, 


o (1) - o (1) 
(2) Pjj Pkk 

Pik ik ~ U — - 

J ’ J S; ~ S k 


(4.67) 


which we rederive below using Wick’s theorem [see Eq. (4.109)]. We can also write 

C 2 } 

Pjkjk i n terms of derivatives of the lRDMs, 

(2) _ ( 2 n ) 2 d 2 T _ 2 n Q Q - (1) 

Pjk, jk — x f)o f)o. — ~T~ f)o \~ r kk/ rjj 


T dsj ds k T dsj 


M) = P?P% + 2n' 


d Sj 


(4.68) 


( 2 ) 

which can be used to evaluate p) k , ]k when Sj = s k or even when j = k. Once we know 
Pp! jk f° r all k 7^ ji an alternative way to find pfj jj is by using the marginalization 
condition 


Y,p { il* = ( 2n - 1 )p. 


(i) 

n 


(4.69) 


Our conclusion is that- the diagonal elements of the 2RDM can be calculated from 
the diagonal elements of the 1RDM. 

The relation (4.67) can be generalized to higher orders; for example, we have 
q i q / 9 


n 

d= 1 


E 


T — Si, \ T S,-, *-.Si,-S jd , 


d= 1 


n 

d'= 1 

d'^d 


(4.70) 


3d 3d 


where j d E {1,2,..., i/}, and j d ^ j d , for d ^ d'. For the general case where the same 
j appears multiple times, a similar procedure can be carried out if required, but it 
becomes increasingly complicated as q increases. 
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In view of the results in this section, we should address the question of whether 
the IRDM encodes all the information about a PCS, including the coefficients A j 
and the Schmidt orbitals ijjj. The answer to this question is a decisive no. One way 
to see this is to recall that the IRDM is diagonal in the Schmidt basis and thus 
is insensitive to the phase of the orbitals. In other words, the IRDM remains the 
same under the transformation 'ipj —> e : ‘'ijjj ; in contrast, the 2RDM and higher-order 
RDMs are sensitive to this phase change. 


4.7 Examples 


In the following, I give some exactly solvable examples which include the case v = 2, 
the totally degenerate case, and the case where the sqs come in pairs. 


For the v — 2 case, we notice that 


(G„ * G„)(x) = 


fX e s i(x~x) e s 2 X 


7= d X 


o Vx'(x-x') 

S 1 + S 2 \ T f 


= 7rexp 


s i — s 2 


X] i 


where 


h(x) = ~ e 

vr ./ 


~ xcos6 dO = — 


x pi —2xu 


: du 


(4.71a) 

(4.71b) 

(4.72) 


7T Jo \/u{l - u) 

is the zeroth-order modified Bessel function.' Putting Eq. (4.71b) into Eq. (4.59), we 
have 


T(s 1 ,s 2 ) = e s+ / 0 (s_) , 


(4.73) 


where 


Si ± s 2 


(4.74) 


Notice that if we added 5 = —(s 1 + s 2 )/2 = —s + to both s x and s 2 , as in Eq. (4.60), 
we would remove s + from T(s!,s 2 ). 

7 

In this dissertation, Ij stands for the jth order modified Bessel function. 
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It is now straightforward to calculate the 1RDM using Eq. (4.47c), 


m 2 n d , , T , ,x 

" rqy 37 , (e +/ “ (s -» 


T 

— n I 1 + 


yxy e s + (ioOO + MO) 

h(s-) 


Similarly, we have 


4 1 ) 


p 22 - nil- 


J o(s-) 


h(s.) 

h (s_) 


(4.75a) 

(4.75b) 

(4.75c) 


(4.76) 


These equations can be compared with fully numerical results, and as shown in 
Fig. 4.3, the two conform quite well in the large-iV limit. 



Figure 4.3: Eigenvalues of the 1RDM (normalized to one) as a function of the 
number of particles. The coefficients X 1 and A 2 are fixed here; i.e., the parameter 
s_ grows linearly in N. The validity of our approximation in the large -N limit is 
confirmed by the numerical results. 


To see the particle-particle correlations in the v — 2 PCS, we calculate the 2RDM; 
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by putting Eq. (4.73) into Eq. (4.47c), we have 

2 d 2 ( s + T , ^ 

* (! - )) 

T(f)s;( e ‘ + ('»(*-)+am)) 


( 2 ) 4n d 
Pii,ii - T j 

2 n 2 5 


= n 2 f _ + 2 Jl ^~^ + i /2 ^~^ 

\2 + Io(s-) + 2 I 0 (s_) 

Similarly, we have 


(2) 2^1 1 -^2 (®—) 

Pi 2 ,i 2 - n ( - - - 


( 2 ) 2 
P'22,22 ~ U 


, 2 I 0 (s_) 

3 Ji(s-) 1 ^(3-) 
2 /o(a-) 2 7q( s -) 


(4.77a) 

(4.77b) 

(4.77c) 

(4.78) 

(4.79) 


It is straightforward to check that these 2RDM elements marginalize correctly, i.e., 
di u 1 + Pi 2 a 2 = Pn and Pigi 2 "h p 22,22 = P 22 The relation (4.67) can be verified 
by using the recurrence relations of the Bessel functions, 


( 2 ) _ / 0 (g_) - I 2 (s_) _ h(s_) _ p ( y L ) - p { S 

P 1 -' 12 2 Iq( s ~) s - Iq( s -) s i ~ s 2 


(4.80) 


In the totally degenerate case (s 1 — s 2 — ■■■ — s v ), all of the eigenvalues of the 
1RDM are the same and thus are determined by the normalization condition, 

Pp = — , for j = 1, 2,..., v . (4.81) 

Similarly, we have that the 2RDM matrix elements pj 333 are the same for all j and 
the matrix elements p-f! ]k are the same for all j 7 ^ k ; moreover, putting the Laplace 
form (4.65) into Eq. (4.47c) implies that 

iorj^k. (4.82) 


We can now use the marginalization condition (4.64c) to determine that 


P { jli k ~ v ( v + 2 ) ( 2S i k + ^ • 


4n 


(4.83) 


We neglect the difference of 2 n and 2n — 1 in the large -N limit. 
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This last result is a special case of the general result for diagonal matrix elements in 
the totally degenerate case. In this situation, plugging Eq. (4.54) into Eq. (4.47c) 
gives us the diagonal matrix elements of the RDMs of all orders, 

(2 n) q d q T(s) 


P 


(<?) 

3l'"3qi 3l"'3q 


T(0) ds 

{2n) q Y{y/2) 


■ • • ds 

3l Ub 3q 


S = 0 


7r 


v/2 


B(qi + 1 / 2 ,... ,q u + 1 / 2 ) 


(4.84a) 

(4.84b) 


n q Y{y/2) 

1 \q + W) 


lh 2 n - !) n 

3 = 1 


(2n) q (u — 2 )!! 
[y + 2q — 2)\\ 


n< 2 % -!)". 

3 =1 


(4.84c) 

(4.84d) 


where q 3 = Yll=i is the number of times j appears in the list of single-particle 

states, j 1; ... ,j q . 


The degenerate case illustrates that when the system possesses special symmetries, 
the expressions for matrix elements can sometimes be solved exactly. Another example 
occurs when the s 3 s come in degenerate pairs, i.e., s 3 = s J+] for odd j. Note that 
Eq. (4.71b) gives the following for = s 2 : 


(G s * G s ) (x) = ne xs . (4.85) 

Using Eq. (4.59) and (4.85), we have 9 

T {xs ) = X~ v/2 ( exp (xsj * exp (ys 3 ) * • ■ ■ * exp (ys^)) (4.86a) 

= y- /2 £ U‘> n -4-) , (4.86b) 

j G odd ^ k G odd 3 ^ 

k+3 

where odd = (1, 3,. .., v — 1}. Note that Eq. (4.86) corresponds to the convolution 
of many exponential distributions, and results from probability theory can be used. 

9 Another way of deriving Eq. (4.86) is by using the residue theorem. 
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The purposes of this section are the following: (i) to find an effective way of approx¬ 
imating the fRDMs; (ii) to gain physical intuition about the fRDMs; and (iii) to 
verify and extend results already obtained. Using Wick’s theorem, 1 will derive a 
relation between the lRDMs of PCSs of 2 n and 2 n — 2 particles. In the large-iV limit, 
this iterative relation can then be turned into a differential equation, and a series 
solution is given for that equation. 

In the Schmidt basis, we have 


— ( vac A n Ojttf, (^) n vac ) 

(4.87a) 

—— Xj\ k ( vac | A n ~ l Gqaj. (^) n 1 vac ) 

n 

(4.87b) 

jr X 3 X k ( vac \A n ~ l (aJ,oq + S jk ) (-4 1 ) " 1 vac ) 

(4.87c) 

4n 2 N n _! / (1) 2 \ 

M (AjAfc (n !) + A jS jk J. 

IM n \ / 

(4.87d) 


This relation implies that p^\n) is diagonalized in the Schmidt basis provided that 
p^\n — 1) is diagonalized. Since p^(l) is diagonalized, mathematical induction 
allows us to conclude that p^\n) is diagonalized. 

That p^in) is diagonalized in the Schmidt basis can also be seen directly from 
Eq. (4.87a): (*4 f ) n | vac ) is a superposition of Fock states that have an even number of 
particles in each of the Schmidt single-particle states, so in the Fock-state superposition 
for a k (*4.^)"| vac ), the single-particle state k always has an odd number of particles; 
thus a k (*A t )"| vac ) is orthogonal to cq vac ) unless j = k. 

For the diagonal elements p'j (//) = Qj[n ), we have 
4n 2 N 

6j{n) = ' n_1 A 2 (p^n - 1) + l) . (4.88) 

Using the condition ]T7 g 3 (n) = 2 n, we can do the iteration without knowing the 
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ratio of the normalization factor N n _ 1 /N n , 


Qj{n) = 2 n 


\) [qj (n - 1 ) + 1 ] 

Ylk= i ^ l [Qk( n ~ 1) + 1] 


(4.89) 


For sufficiently large n, all the probability concentrates on the dominant eigenvalue 
A]. To get useful results, we again assume the differences of the A^s are small and 
use the parametrization (4.43), A 2 = 1 + Sj/n. To turn the iterative equation into a 
continuous differential equation, we introduce the parameter r E {1/n, 2 jn ,... , 1} 
to represent the iteration steps. In addition, we defer normalization to the end of the 
process, instead of imposing it at each iteration. The new iterative equation reads 


8 1 ^ 

Qj(nr + 1) = (l + — ) ( 'Qj{nT ) + — ^ £fc( nr )) (4.90a) 

T k =1 

1 i v 

- Qj(nr) + - ^ SjQj^nr) + — ^2 Qk{ nr ) ) , (4.90b) 

k =1 

where the term of order 1/n 2 is neglected, and we use (l/2nr) g k instead of 1 to 
deal with the fact that o 3 is unnormalized. Taking the limit n —> oo, we have the 
differential equation 



k= 1 


(4.91) 


The problem with Eq. (4.91) is that it diverges at small r due to the factor l/2r. 
This divergence is an artificial consequence of our decision to defer the normalization 
to the end; one remedy is to modify the differential equation to 



fc=i 


(4.92) 


where the extra term, which only introduces an overall factor, keeps g 3 from diverging 
from our initial condition £?j( 0 ) = 1 , for j — 1,2,... ,u. 


For the case v = 2, we can decouple the two equations in Eq. (4.92) by introducing 
g± = Qi ± q 2 ] the modified Bessel functions / 0 (t) and /i(r) solve these equations. 
Thus we have recovered our former results (4.75) and (4.76) using an entirely different 
approach. 
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For the general case v > 2, the solution of Eqs. (4.92), which are linear equations, 
can be expanded into a series. This can be done most conveniently by introducing 
the matrix equation, 

^-=( z ~^(Z-S))nT), ( 4 . 93 ) 

where T(r) is the transition matrix, Sj k — 1/u (for j,k — 1 , 2 ,,u) is the projector 
onto the symmetric vector (1,1,..., 1) T ; Z = diag(s,, s 2 , ■ • •, s v ), and 1 is the identity 
matrix. Knowing the transition matrix allows us to solve for 

V V 

Qj(t) = Y T jk ( j) Qk{ 0) = Y T A r ) > ( 4 - 94 ) 

k=1 k=1 

where the initial condition is p fc (0) = 1 (for k = 1, 2 ,,u). By choosing the initial 
condition 10 of the transition matrix to be T(0) = S, we have 

OO 

T(r) = ^+^r m T m , (4.95) 

m =1 

where T m , m = 1,2,..., are matrices to be determined. Putting Eq. (4.95) into 
Eq. (4.93), we get 

OO OO OO 

Y, m T ™~ 1 T m = Z S +Y Tm ZTm~\Y Tm_1 ( 1 - S ) ■ ( 4 ‘ 96 ) 

m =1 m =1 m= 1 

Comparing the coefficients of different orders of r in Eq. (4.96), we have 


(m 1 + ^ (l - 5)) T m = ZT m _ x , for m > 1 . 
Because the matrix m 1 + u(l ** <S')/2 is invertible, we have 


m _l + ^/ 2 m vrr 

J- m A -*■ m— 1 


(4.97) 


(4.98) 


m + v/2 

Note that the solution to Eq. (4.98) always take the form T m = Poly (Z) S, where 
Poly(Z) is some polynomial of the Z matrix . 11 This is because the identity 


uSDS = tr (D) S (4.99) 

10 T 

We can choose any matrix that stabilizes the vector (1, 1,..., 1) to be T(0). The 

choice T(0) = S is most convenient, because the term of order r 1 in Eq. (4.96) disappears 

automatically. 

11 For example, we have T 0 = S. 7j = 1+ ^ 2 S , and T 2 = (i^_ t )/ 2 )( 2 +' > // 2 ) w here we 
assume tr(Z) = 0 without loss of generality. 
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holds for any diagonal matrix D. Thus the matrix equation (4.98) can be turned into 
a c-number equation for polynomials. In addition, we see that the matrix norm of 
T m begin to fall quickly after m > s 1; which gives us an estimate of how many terms 
we need to include in the expansion to get a desired precision. 

In this section, we have discussed an alternative way of deriving the lRDMs with 
a differential equation approach. The results we have here conform with those derived 
previously. The solution to the differential equations can be expressed as an infinite 
sum. Compared to our other approaches, this one allows one us to solve the lRDMs 
efficiently, but approximately, and is suitable for numerical evaluation. 


4.9 Determining the 2RDMs 

In the former sections, we only discussed the diagonal elements of all RDMs. Here, 
we show how to represent the off-diagonal elements of the 2RDMs with their diagonal 
elements. In the Schmidt basis, the 2RDM reads 

pS 2 , jlj2 = ^( yac \ A ™ 44^2^1 (^) n I vac > • ( 4 - l0 °) 

All the matrix elements of p ^ are real, and it satisfies the normalization condition 

tr p< 2 > = £ = 2n(2 „ _ j) . (4,101) 

3,k 

Using Wick’s theorem, we find the 2RDM must have the form 

^ hk^fak-^ j (4.102) 

where the real matrices £ and £ 7 are symmetric to ensure the Hermiticity of the 
2RDM. 

In the large -N limit, this form can be further simplified by the condition Xj/Xk = 
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1 + 0(1/n) [see Eq. (4.43)], 


N (vac M a^a^a^al^A^) vac) 

(4.103a) 

N (vac \A al^a^a^ (A*) vac ) + 0(n) 

(4.103b) 

Kl ~A^ VaC a C a h ak 2 a h (^ f ) n vac ) + 0(n) 

(4.103c) 

^ ( vac \A n a\ i a) 2 a k 2 a ji (M T )” vac ) + 0(n) 

(4.103d) 

Phk 2 ,k lj2 + 0(n) , 

(4.103e) 


f 2) • • 2 

which means p\‘ k j]]2 is symmetric at the leader order, n , under the exchange of j 1 

and k\. Thus we can further simplify Eq. (4.102) into 

Ppk 2 ,j 1 j 2 — ^■jlk 1 ^jd2^ k l k 2 + (^'l fc l^i2 fc 2 + ^3l k 2^h k l ) ’ (4.104) 

where the matrix £ determines p^. Note that all the elements of £ are nonnegative 
due to the positiveness of the Schmidt coefficients. Moreover, the positiveness of p ^ 
is equivalent to the positiveness of the matrix £ jfc + 2h jfc £ jj -. Conversely, the matrix £ 
can be determined by the diagonal elements of p^ 2 \ 

(4.105) 

and these diagonal elements can be determined from Eq. (4.47c). 

In Sec. 4.6, we derived a relation (4.67) between the diagonal elements of p ® and 
those of p^ l \ Here we rederive this relation using Wick’s theorem, without making the 
large -N approximation. By doing contractions of the single annihilation and creation 
operators with the pair creation and annihilation operators, we have for j ^ k, 


N pfkjk = (vac a ] a l a k a j (^ t ) n |vac ) (4.106a) 

= 4n 2 A 2 ( vac \ A n ~ 1 aja\a k a}^ (^) n *| vac ) (4.106b) 

= 4n 2 A 2 ( vac \A n ~ 1 (a^a\a k aj + a\a k ) (A^) n X | vac ) . (4.106c) 

Another way of contracting gives us 

N pf k j k = 4n 2 A 2 k ( vac | A n ~ 2 a k a^aja\ (A^) n 1 \ vac ) (4.107a) 

= 4n 2 A% ( vac \ A n ~ 1 (ajaj,a fc a J - + (^) n X | vac ) . (4.107b) 
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Multiplying Eqs. (4.107b) and (4.106c) by A^ and X 2 k . respectively, and subtracting 
the results gives 

(A 2 — A l) N pflj k = 4n 2 \ 2 \ 2 k (vac | A n ^ x (a]a_j — a\a k ) (.A 1 ")™ 1 | vac ) (4.108a) 

= ( vac |*4“ (A l Ojdj — \ 2 a\a k ) (^)"| vac ) , (4.108b) 


which leads to 


,2 ( 1 ) \ 2 ( 1 ) 
(2) _ A k Pjj A j p kk 

Pjk,jk ~ 


JJ _ 

A?-A l 


(4.109) 


With the large -N parameterization (4.43), Eq. (4.109) reproduces Eq. (4.67) in the 

( 2 ) 

large -N limit. This relation allows us to determine p jk for j ^ k, from the 1RDM, 
( 2 ) 

and then p 33 r] is determined by the marginal condition, 


ppjj = ( 2n - 1 )pp 


(4.110) 


In the large -N limit, we have determined the 2RDM by using the 1RDA1 in the 
Schmidt basis, which suggests that the 2RDM is a function of the 1RDM. This 
conclusion assumes, however, that the Schmidt basis can be calculated from the 
1RDM. To be explicit, given the 1RDM, one can diagonalize it and thus find the 
Schmidt coefficients and Schmidt orbitals up to an arbitrary phase for each orbital. 
The 2RDM (4.104), written out explicitly in the Schmidt basis, takes the form 

P {2) = ^^jk(\^k){^j I ® \^k)(^j I (4.111a) 

j,k 

+ + ■ (4.111b) 

The terms in Eq. (4.111a) are sensitive to the phases of the Schmidt orbitals when 
j k and thus cannot be determined from the 1RDM. 

In the case v = 2, we can solve the 2RDM exactly in the large-fV limit. Using 
the results (4.77), (4.78), and (4.79) and also Eq. (4.104), we have the following 
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expression for the 2RDM in the Schmidt basis, { 1 11 ), 1 12 ), | 21), | 22) }, 
3/ 0 + 411 + J 2 0 0 I 0 — / 2 



0 Iq — 12 -^0 ~~ h 0 
\ Iq ~ h 0 0 3J 0 — 4 I x + I 2 J 


0 



0 


(4.112) 


where J 0 , /,, and J 2 are the zeroth, first, and second order modified Bessel functions 


with argument s_ = — s 2 )/2. Equivalently, we can write in the Pauli basis 



(4.113) 


The two-particle state (4.112) is not entangled; i.e., it has zero concurrence. It is 
known that all pair-wise entanglement vanishes in large bosonic systems due to the 
monogamy of entanglement [KBIOO]. 

Another case that can be solved analytically is the totally degenerate case where 
Si — s 2 — ■ ■ ■ — s v . Using Eqs. (4.83) and (4.104), we have 



(4.114) 







Chapter 5 


Applying Pair-Correlated States to 
Fragmented Condensates 


The fundamental laws necessary for the mathematical treatment of a large 
part of physics and the whole of chemistry are thus completely known, 
and the difficulty lies only in the fact that application of these laws leads 
to equations that are too complex to be solved. 

- Paul Dirac 


In the last chapter, we introduced the Bosonic Particle-Correlated State (BPCS), a 
state of N — l x n bosons that is derived by symmetrizing the n-fold tensor product 
of an arbitrary /-boson (pure or mixed) state <r®. In particular, we are interested in 
the pure-state case for l = 2, i.e., cr (2) = | |, which we call a Pair-Correlated 

State (PCS). Note that PCS can also be generated by projecting a multimode squeezed 
vacuum to the iV-particle sector . 1 When there is one dominant squeezing parameter— 
or equivalently, a dominant Schmidt coefficient in the two-particle wavefunction—PCS 
reproduces the number-conserving Bogoliubov approximation. For the case where 
many squeezing parameters are of the same size, PCS describes a fragmented state 

1 With the “quasi-normal-ordered” factored form of the squeeze operator [Per77, Hol79], 

see Eq. (4.28), the squeezed vacuum takes the form | !^ sv ) ~ exp (—4 tanhy fc aj, 2 ) | vac). 
Projecting this state to the N particle sector gives the PCS. 
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with nearly maximized two-particle quantum correlations. 

In this chapter, I will discuss how to use the results from the last chapter to derive 
equations that determine the PCS ground state. I also formulate time-dependent 
equations for the evolution of PCS using the technique described in Chap. 2. For the 
two-site Bose-Hubbard model, we calculate the trace distance between the 2RDMs 
given by the full numerical solution and the closest PCS. The results confirm that 
PCS provides a good representation of the ground state over the entire parameter 
space; more interestingly, for the time-dependent case our numerical simulations 
suggest that the error in the 2RDM given by PCS does not get larger with increasing 
evolution time. 


5.1 PCS Ground State 

In this section, we discuss how to solve for the PCS ground state using a variational 
method. By fixing the PCS parameters s (the eigenvalues of the 1RDM are thus 
fixed), we find equations that determine the orbitals. By fixing the orbitals, we write 
the energy as a function of s, and we find the s that minimizes the PCS energy. 
Repeating this whole procedure several times allows one to find the PCS ground state 
approximately. We also derive a condition for fragmentation using the PCS ansatz. 

Consider the state of N bosons modeled by the following Hamiltonian, 

H = 'H X + H 2 

= j ^ t (x)(-^-V 2 + R(x)jrJ;(x)c/x+| J [rjj t (x )] 2 r|) 2 (x) dx , 

where FL X and FL 2 are the single- and two-particle Hamiltonians respectively. The 
energy expectation value of the Hamiltonian (5.1) for a PCS reads 

£pcs = (^pcsl^l ^pcs ) = < ^pcs I (ft 1 + n 2 ) I %cs > • (5.2) 

The contribution from the single-particle Hamiltonian is 

(^pcsl^il^pcs) =tr(p (1) /7 1 ) = ^Qj / V'K X )(-^—V 2 + R(x))^( x )dx , (5.3) 

A — -\J — M 
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where o 3 = is the occupation number of the jth orbital ^y(x); the contribution 
from the two-particle Hamiltonian is 


(^pcsl^pcs) 


| J ( ^pcs I [^ t (x)] 2 iJj 2 (x) | F pcs ) dx (5.4a) 

£ v n (2) 

2 / v rlm,jk 

j,k,l,m=l 


^i(x)^fe(x) ^ m (x)^,(x) dx 


(5.4b) 



where we used the relation (4.104). Note that the overall phases of the Schmidt 
orbitals are important in the first term of Eq. (5.4c). Combining Eqs. (5.3) and (5.4), 
we have 


Epcs ^ Qj 
3=1 


in 


R , 


2m 


€ +v + ^m + tti Qi^Wi + Qmt 


' 3 r jr j 


dx , (5.5) 


where the effective potentials are given by 


#i( x ) = — ( 2 + s jk)£jk l^fc( x )l J , <5i( x ) = — ^( x ) • ( 5 - 6 ) 

e 3 k =1 6 * k= 1, 

Ml 

Note that A and Q are of order IV 0 , because p is of order IV -1 , is of order IV 2 , 
and Qj is of order IV. 


The energy expectation value A pcs depends both on the Schmidt orbitals and the 
PCS parameters s. As a first step, we only allow the Schmidt orbitals to vary by 
fixing s. To keep the orbitals orthonormal, we introduce the Lagrange multiplier 
L = k dj k ('ijjj 7/4 ), where 'd ]k = x)* kr Varying the jth orbital, we have the following 
equation for the PCS ground state, 

_ 1 $(.E pcs - L ) 

Qj H*j ( x ) 

/ ih 2 \ u 

= (-^v 2 + ^( x ) + ^( X )J^( X ) + Qj( x )^( x ) - 5^tfjfcV’k( x ) • ( 5 - 7 ) 


Combining the above equations for j = 1,2,... , 1 /, we can determine the optimal 
orbitals; a promising candidate to solve these differential equations numerically is the 
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imaginary-time evolution method, where at each time step i) is chosen so that the 
evolved orbitals remain orthonormal. 


After the orbitals are fixed, the energy of PCS takes the simple form 

v i V 

E pcs = 6jQj + - (1 + 28jk)£jkrijk ) ( 5 - 8 ) 

j= 1 j,k =i 

where 

e j = I ^*( x )(-^ v2 + l/( x ))^( x ) dx , (5.9) 

V j k = Y+25~ k J ([^( x )] 2, 0fc( x ) + 2 |Vh( x )| 2 |^ fc ( x )| 2 ) dx • (5.10) 

Note that the elements of rj ]k are determined by the Schmidt orbitals and generally are 
sensitive, through the first term in the integral, to the overall phases of the Schmidt 
orbitals. For the PCS ground state, the chemical potentials of different orbitals must 
be the same 


dE pcs 

d Qj 


e j+ 9 ^ + 25 tk) 


j,k =1 


dijk , t drj jk \ 
doj Vjk + ^ jk d L>j ) 


fl j 


where the partial derivative d/dQj is taken by fixing Q k for all k ^ j. 


(5.11) 


Here we discuss the case v = 2 in detail. In this case we can assume that ?/ 12 
is real for the following reason: with everything else held fixed, we can lower the 
energy of PCS by varying the overall phases of the two orbitals so that the integral 
J[*01 ( x )] 2 *02 ( x ) dx is negative, thus making r/ 12 always real for the PCS ground state 
and equal to 

V 12 =9 j |Vh(x)| 2 |^ 2 (x)| 2 dx . (5.12) 

Note that this argument does not work for higher v because there are several com¬ 
peting rj jk s to consider. Moreover, for u — 2, one can argue that the phases of 
the two orbitals do not depend on x, because the kinetic energy term and the in¬ 
tegral J[?/q(x)] 2 0 2 (x) dx can be minimized simultaneously by making the phases 
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homogeneous. 2 Using the Cauchy-Schwarz inequality, we have 

V 12 =9 J |Vh(x>| 2 |^ 2 (x)| 2 hx < V 7 n 722 ■ (5.13) 

With the results (4.75) and (4.76) from Chap. 4, we have 

Qi = n (1 + Xi) , Q 2 = n(l-Xi) , (5-14) 

where yq(s_) = Ii(s_)/I 0 (s_) G [0, 1) is the ratio of the first and zeroth order 
modified Bessel functions, with s_ = (.sq — s 2 ) /2 being the single relevant PCS 
parameter. Without loss of generality, we assume that the occupation number 
of the first mode is greater or equal to that of the second, i.e., s_ > 0. Using 

Eqs. (4.75)-(4.79) and the relation (4.105), we have 

2 2 2 
77 / 77 / 77 / 

£n = (3 + 4 Xi + X 2 ) 1 £12 = (l — w) ’ ^22 = (3 4xi + X 2 ) , (5.15) 

where X 2 ( s -) — 7 2 (s_)// 0 (s_) G [0, 1). Note that the elements of as functions of 

s_, are properties of how the Schmidt orbitals are populated within the PCS. 

The energy expectation value thus reads 

£pcs ^1 Qi + e 2^2 + 2 + 3^22 722 + £ 12(712 + 712 ) j (5.16a) 

= E 0 + (n(e 1 - e 2 ) + n 2 (r; n - j] 22 ) Jxi + ~j{vu + V 22 ~ ^ 12 )X 2 (5.16b) 

= E 0 + n(-CiXi(s_) + 77 X 2 ( 0 ) • (5.16c) 

Here 

n 2 

E 0 = n (e x + e 2 ) + -^-(377 n + 3rj 22 + 2r] 12 ) (5.17) 

is the energy expectation value for equal occupation numbers, g 1 = q 2 (i.e., s_ = 0); 

ci = (e 2 + np 22 ) - (ei + n^n) (5.18) 

is the difference in chemical potentials at equal occupation numbers (with the orbital 
wavefunction fixed); and 

77 / 

C 2 = 2 fall + d 22 - 2712 ) (5.19) 

2 

This also might not be true for v > 2, since the minimization of the integral 
/[Vh (x)] 2 i/j 2 (x) dx can become x-dependent and thus compete with the kinetic term. 
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measures the difference between the strengths of intra- and inter-orbital interactions. 
From the condition (5.13), we have that c 2 > 0. Note that in the second term in 
Eq. (5.16c), rq and c 2 contain the dependence on the Schmidt orbitals, and Xi and 
Xi contain the dependence on how the orbitals are populated within the PCS. 

Since dE pcs /ds_ = n (dE pcs /dg x — dE pcs /dg 2 ) dxi/ds_, the condition (5.11) that 
the two orbitals have the same chemical potential is simply the extremal condition, 
dE pcs /ds_ = 0, satisfied by the minimal energy. Since y^s.) is an odd function and 
y 2 (s_) is an even function, the minimum of E pcs (s _) occurs when s_ has the same 
sign as cq. If cq is negative, the minimum has s_ < 0, violating our assumption that 
Qi A f? 2 i this is telling us that we should switch the roles of the two Schmidt orbitals 
to make cq nonnegative. In what follows, we assume that cq > 0 consistent with our 
assumption that Qi > g 2 . 

The function E pcs (s _) seems to be monotonically decreasing when cq > 2c 2 , 3 so 
the PCS ground state is a pure condensate, i.e., s_ = +oo. When cq < 2c 2 , the terms 
with coefficients cq and c 2 in Eq. (5.16c) compete with each other, and the PCS ground 
state is a fragmented condensate. In Fig. 5.1, the occupation ratio Q\/N for the PCS 
ground state is plotted as a function of cq/c 2 ; the data are calculated numerically. 
From the figure, we see that the particles are distributed evenly between the two 
orbitals when cq/c 2 = 0, and there is a single condensate when cq/c 2 = 2. Note 
that this single-condensate result does not conflict with the claim made in Sec. 4.4 
that PCS reproduces the number-conserving Bogoliubov approximation, because the 
depletion of a single condensate predicted by the Bogoliubov approximation is of 
order 1/N and goes away in the large -N limit. Given these considerations, we have 
the following condition for fragmentation, 

Ci — 2c 2 = e 2 — e 1 + N(r) 12 — r) u ) < 0 ; (5.20) 

fragmentation happens when the two orbitals are close to degenerate, and their 
wavefunctions have little overlap. 

3 

For large s_ we can prove this rigorously by using the following asymptotic expansion 

of the modified Bessel functions: I Ax) ~ %— ( 1 — ^ + (4q ~ 9 ^ + ■ ■ ■ Y 

; P2F x \ 8x 2!(8x) J 
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Figure 5.1: The occupation ratio Qi/N for the PCS ground state plotted as a 
function of cq/c^. We have equally populated orbitals for c x jc 2 = 0 and a pure 
condensate for c x fc 2 = 2. 

If we choose the dominant orbital t\) x to be the GP ground state </>, then there is 
no -02 such that the condition (5.20) can be satisfied for a single-component BEC 
in the large -N limit. This can be proved by noticing that the quantity in Eq. (D.5) 
is nonnegative; by choosing i/j 1 (x) = </>(x) to be real and yj 2 (x) = </q_(x) to be pure 
imaginary for all x, one has that e 2 — e x + N(r] 12 — r] u ) > 0. Since 0j_ can be any 
single-particle state orthogonal to the condensate wavefunction, the condition (5.20) 
cannot be satisfied. The dominant orbital of the PCS ground state, however, does 
not have to be the GP ground state, and a different choice could lead to a fragmented 
BEC. 

The two intertwined aspects of describing a many-body system such as a BEC 
are, first, which single-particle states to populate and, second, how to populate them. 
Within the PCS formalism, the choice of single-particle states is taken up by the 
choice of Schmidt orbitals; PCS makes a particular choice for how to populate these 
orbitals, with the remaining freedom within the formalism being the choice of the 
Schmidt coefficients. For the case v = 2, the dependence on orbitals is contained 
in the quantities iju, 7] 22 , and ?/ 12 , and the freedom in populating the orbitals is 
described by the single parameter s_. We can study the effectiveness of the PCS 
ansatz as a way of populating single-particle states by considering contexts in which 
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all the relevant properties of the single-particle states are taken into account by just 
a few parameters appearing in a many-body Hamiltonian. 

One such context, which we turn to now, is a two-site Bose-Hubbard model with 
Hamiltonian 

^tbh — — r(a\a 2 + (4 a i) + 2 (A/'ifVV’i — 1) + J\f 2 (J\f 2 — 1)j , (5-21) 

where r is the tunneling strength, A/} = aja,,- is the particle number operator of 
the j t-h mode, and u > 0 describes the on-site repulsion. We numerically calculate 
the exact ground-state vector | F exc (u/r)) as a function of u/r for N = 200 bosons 
using the imaginary-time evolution method; we then calculate the exact single- and 
two-particle RDMs as a function of u/r using the ground-state vector. We 
construct a PCS 2RDM p^} s from Eq. (4.112) by assuming that its marginal is the 
exact 1RDM pX ) c : the relative phase of the two Schmidt orbitals, which p^X does not 
provide, are chosen to minimize the energy. 

The error of the PCS 2RDM, measured by the trace distance \ tr|pj/2, — pixel) is 
plotted in Fig. 5.2 as a function of the occupation ratio Pi/N (which depends only on 
u/r); the error is smaller than two percent over the whole parameter space (hardly 
any error would be seen had we used fidelity instead of trace distance). We compare 
the PCS result to a widely used ansatz, the double-Fock state (DFS), where and 
N 2 particles occupy the two orbitals. We also plot, in Fig. 5.2, the trace distance 
§ t r |Pdfe — Pexcl minimized over all p^ that are consistent with the DFS ansatz, as a 
function of the ratio Qi/N. This comparison shows that PCS provides a much more 
faithful representation of the ground state of the 2-site Bose-Hubbard model than 
DFS. 

Before diving into the time-dependent equations for the PCS ansatz, let us 
look at how faithfully PCS can represent the time evolution of the two-site Bose- 
Hubbard model. We consider the time evolution of the 2-site Bose-Hubbard model of 
N = 1000 atoms. Initially, all the atoms are in the ground state of the noninteracting 
Hamiltonian (u = 0, r / 0), and then strong on-site interactions are suddenly turned 
on (i.e., u is made very large compared to r). We calculate the evolution of the whole 
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Figure 5.2: For the ground state of the 2-site Bose-Hubbard model of iV = 200 
atoms, we plot the trace distances (normalized to one) of the exact 2RDM to the 
closest 2RDMs given by PCS and DFS as functions of the occupation ratio Q\/N. 
The particles are split into two halves when Q\/N = 0.5 {u/r —> oo), whereas all 
the particles are condensed into a single mode when Qi/N = 1 (u/t 0). The 
small error at Qi/N = 0 is due to the finite number of particles and goes away for 
N —» oo. 


state vector numerically using the Crank-Nicolson method, which has second-order 
precision and is always numerically stable; the exact 2RDM is derived from the state 
vector. We then optimize the parameter spaces of PCS, DFS, and a Gross-Pitaevskii 
State (GPS) to minimize their trace distances to the exact 2RDM. The errors of the 
2RDMs, measured by the trace distances normalized to one, are plotted as functions of 
evolution time in Fig. 5.3. The oscillation of the error given by GPS is a consequence 
of the collapse and revival of phase [GMHB02, SDZ11]; i.e., the purity of the 1RDM 
oscillates. These numerical results suggest that PCS might be useful to describing 
the dynamics in the strongly interacting regime. 


5.2 Time-Dependent Equations for PCS 


In this section, I discuss how to derive time-dependent equations for PCS using the 
“projection” technique developed in Chap. 2. The key idea is to evolve the 1RDM 
under the full Hamiltonian and then update the PCS using the evolved 1RDM. This 
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Figure 5.3: We consider the time evolution of the 2-site Bose-Hubbard model 
of N = 1000 atoms. Initially, all the atoms are in the ground state of the 
noninteracting Hamiltonian (a = 0, r ^ 0), and then a strong on-site interaction 
is suddenly turned on. We plot the trace distances (normalized to one) of the 
exact 2RDM to the closest 2RDMs given by PCS, DFS, and GPS as functions of 
evolution time. 


method, however, is insensitive to the crucial relative phases of the orbitals, so we 
calculate these phases separately by directly evolving the state vector. The result 
is a set of coupled equations that determines the evolution of the orbitals as well 
as the occupation numbers for these orbitals. Our results, a generalization of the 
time-dependent GPE, enables one to solve for the dynamics of a fragmented BEC. 

Suppose we have a PCS at time t, 

\^pcs(t)) = ^[^Wrivac) , A\t) = [a](t)] 2 , (5.22) 

where aj(t) = (i[> | ipj (t )) is the creation operator of the jth orbital; these orbitals 
^■(x, t), for j = 1,2,, oo, form a complete, orthonormal basis of single-particle 
states (note that we are allowing j > v). We do not put any time dependence into 
N, because such time dependence can be absorbed into A\t) as a rescaling of the 
Schmidt coefficients A j(t). 
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Consider now the state | & pc8 (t )) evolving under the BEC Hamiltonian 

n(t) = n 1 (t) + n 2 (t) ( 5 . 23 a) 

= f 4 >t (x)(-^V 2 + H(x,t) jrKx)dx+^ j [x|d(x)] 2 r^ 2 (x) dx , (5.23b) 

where Fi 1 and 'Ho are the single- and two-particle Hamiltonians respectively. After 
letting the state (5.22) evolve for a time dt under the Hamiltonian (5.23), we “project” 
it back into the PCS submanifold; the rate of error in the state vector that the 
“projecting” step introduces is 

C(t)) = 'Ad I % cs(0 ) - I HW ) , (5.24) 

where | F pcs (t )) stands for the time evolution within the manifold of PCS states. 


As in the product ansatz developed in Chap. 2, minimizing the error introduced 
by the projection means requiring that the error vector | d/ err [t)) be perpendicular 
to any variation within the PCS manifold. Thus the evolution of the PCS state is 
determined by the condition 


(<^pcs (t) I ^err (t )) = 0 for all 5F (t) . 


(5.25) 


For a small variation of the two-particle creation operator A^(e) = A^ + eV^/n, where 
V can be any quadratic function of the field operator i[>(x), we have 


M p CS ) = ^|^p C s(e)) 


e=0 




(^ t ) n “ 1 |vac) . 


(5.26) 


Note that we treat the normalization factor N as constant when A changes; i.e., we 
allow the perturbed state to be unnormalized. Putting Eq. (5.26) into Eq. (5.25) 
gives 

( vac | [^4(t)] n_1 V | 'i'err(f) ) = 0 for all V quadratic in . (5.27) 

Putting Eq. (5.24) into Eq. (5.27) determines the evolution of the PCS state | ^ pcs ). 


In particular, for V = A, we have | SF pcs ) = | !^ pcs ); putting this into Eq. (5.25) 
and keeping only the real part gives 

o = Re (^ ($pcs (t) | H{t) I & pcs (t )) - (tftpcs (t) I ^p CS (t)) ) (5.28a) 


= - Re ((^ pcs (7) | !Ppcs(t))) , 


(5.28b) 
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which says that our procedure preserves the norm of | ^ pcs ), 
(tfpcsWI V pcs (t)) = constant . 


(5.29) 


We can write Eq. (5.27) equivalently as 

( vac | [^4(f)] n_1 a,j(t)a k (t ) | F eir (t) ) = 0 , for all j, k e {1, 2, ..., oo} , (5.30) 


which will be used throughout the rest of this dissertation. Directly evaluating 
I ^pcs(^)) with Eq. (5.30) is hard, because it usually involves calculating the three- 
particle RDM arising from the product of the 'H 2 and the term cij(t)a k (t). To avoid 
this complication, we notice that for all j, k G {1, 2,..., oo}, 


(tfpcsW I a 5(*) I ^errO) > 


-^= Aj-(t) < vac | [M(f)] n 1 a_j(t) a k (t) \ & err (t ) ) = 0 , 

(5.31) 


where = 0 for all j > v. This equation can be recast in terms of the field operators 
as 

OO 

0 = X^ Vh*(x, t) ^fc( x , t) ( ^pcs (t) I a](t) a k (t) | ) (5.32a) 

j,fc=i 

= < ^pcsW I ^ T (x) rj)(x) | \P evi {t) ) . (5.32b) 


The physical content of Eq. (5.32b) is that the 1RDM, 

P (1) (x, x ; t) = ( ^(t) | tb t (x / ) t|>(x) | ^ pcs (t) ) (5.33a) 

OO 

= ^*( x ’ t) ^fc( x , *) ( ^pcsW | a i(^) <**(*) | ^pcsW > (5.33b) 

3,k =1 

V 

= X] e i ^j( x > 0 ( x > *) ’ (5.33c) 

j=i 

is left unchanged by the projection procedure; i.e., the change in the 1RDM under the 
evolution within the PCS manifold is the same as the change under the exact evolution. 
Thus a promising procedure for determining the dynamics of PCS is to evolve p^(i) 
exactly for a short time dt and then to update | F pcs (t + dt)) using p^\t + dt ). The 
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problem with this procedure is that the 1RDM does not encode the phases of the 
orbitals; i.e., p ^ remains the same for the transformation ^(x) — y e* 6 ' i '0j( x ) with 9j 
real. These phases, however, play a crucial role in determining p^; for example, the 
process of annihilating two particles in the jtli orbital while creating two particles in 
the fcth orbital feels the relative phase e^ ej ~ 9k \ 

Despite this problem, we can make substantial progress by considering the evolu¬ 
tion of the 1RDM. According to the BBGKY hierarchy, p ^ is a function of p when 
there are only up to two-particle interactions; luckily, we anticipated this situation 
and considered the 2RDM in detail in Chap. 4. The exact evolution of the 1RDM 
reads 


iftp (1) (x,x'; t) = (F p „,(*) | [\Jj t (x')\Kx), H(t)\ | ^ pcs (t) ) . (5.34) 

The commutator in Eq. (5.34) can be divided into two parts: 

[ip t (x')^(x), H(t)\ = [-ilb(x')xKx), fti(t)] + [^ t (x / )^(x), n 2 (t)\ . (5.35) 


The first part, due to the single-particle energy, corresponds to a common unitary 
evolution of the orbitals. The second part, due to the interparticle interaction, assumes 
the form 


[V(x')M>( x ), % 2 ] 


9_ 

2 


^ T (x)^(x), [\|k(x")] xp 


t 2 ii, 2 (' x "') 


dx" 


9 (^ t (x)^ t (x)^ 2 (x) - [^ T (x , )] 2 iKx / )iKx)j . 


(5.36a) 

(5.36b) 


Putting Eq. (5.36) into Eq. (5.34), we have the following time derivative of the 1RDM 
induced by Fi 2 , 


ifoPul = ( ^pcs | Wx) tKx), H 2 \ | ^pcs > 




■ pcs 


4 )T (x / ) 4 * T (x) 4 )2 (x) - [il) t (x , )] 2 rKx , )rKx) 




9 p( iZ jk ( X 0 ( X ) V’m (x) - V’fc (x) 1 pm (x 

j,k,l,m =1 


pcs 


X . 


(5.37a) 

(5.37b) 


(5.37c) 
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Using the large -N expression (4.104) for the 2RDM, we have 

V 

ih PH 2 -9j2 ( £>jl ^jk^lm ^sjkiPjl^km ^jm^kl) ^ 

j,k,i,m= i (5.38a) 

X V'l(x') (^fc( x ) l(x) - V’fc(x') V’m( x/ ))v , z(x) 

V 

= g k (^( x ) ^*( x ) ^( x )- V’jW ^j( x ) [^(x )] 2 

(5.38b) 

+ 2|^ fc (x)| J ^(x)^*(x / ) -2^(x)Vh( x ) l^fc( x, )| 2 ) • 

Introducing the same effective potentials as were defined in Eq. (5.6) for the time- 
independent case, 

V V 

Rjt*) = — ( 2 + S jk)tjk l^fc( x )l"’, <2j( x ) = — V’fcW , (5-39) 

6 i k= 1 ^ fc=l, 

kAj 

we have the following time-dependent equation for the evolution of the Schmidt 
orbitals: 

• / ft 2 \ 

ih'l’j = ( - ^ v2 + u + A, J Vb + Qjrf for j G {1,2,..., v} . (5.40) 

Note that Eq. (5.40) reduces to the GP equation when v — 1; for y>lwe have v 
coupled nonlinear equations describing the dynamics of the orbitals. 


As mentioned before, our procedure of using the 1RDM to update the PCS fails 
to capture the dynamics of the phases, Im( Cj \i J j)- Moreover, Eq. (5.40) does not 
maintin the orthonormality of the orbitals. I turn now to discussing how to fix these 
two glitches. 


To evaluate the dynamics of phases of the orbitals, we use the condition (5.30), 


0 = (vac | A n 1 ajd k | F eri .) 


vac 


ATI— 1 

CLj CL 




(5.41) 


which is stronger than the approach of evolving the 1RDM and can be used to 
determine the orbital phases. In the second-quantized picture, a change of the phase 
of the jth orbital takes the form ihoj = fij a], where the real number pj is to be 
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determined. The dynamics of the pair creation operator for this restricted change 
reads 


d 


^ = vX - 2 S x A a ) = ^ vk f > ( 5 - 42 ) 

j =i 3= 1 j=i 

which is diagonal in the Schmidt basis. Equation (5.41) can thus be simplified to 
Im ( vac 


Re ( vac 
d 


d\, 


A n 1 cijCij — (hi (A^) n — /j, k X k a\.a\(A^) n vac^)=0 (5.43a) 

fe=i 

V 

AC -1 fljgj (% — Pk a \ a k') (^) n vac ^> = 0 (5.43b) 

fc=i 

.A n hk a k a k) (^)” vac ) = 0 • (5.43c) 


vac 


k=1 


Thus the expectation value in Eq. (5.43c) remains unchanged for any infinitesimal 
increment dX of the coefficients A. When dX is parallel to A, this expectation value is 
simply rescaled, which gives the condition 


vac 


A n (n~ ^ /ifc4 a fc) (- 41 ) 


k=1 


vac ) = 0 


(5.44) 


For any dX keeping the normalization factor Nr unchanged, we have 



Ft ^ ^ 


k= 1 


^pcs(A) ) = d (e( A) - dkQk( x )) 

k=1 


0 . 


(5.45) 


Because we have the constraint J2k =l ^Qk — 0, this result implies that p k equals the 
chemical potential of the kth orbital up to a constant, which can be determined by 
using Eq. (5.44). The constant, however, only adds an overall phase to the PCS state, 
and we can neglect it by letting pj be the chemical potential of the jth orbital, 


Pj 


0E(q) 

d Qj 


fixing all Q k s for k A j , 


(5.46) 


where E(q) — E(n, A) is the energy expectation value. According to Eq. (5.40), the 
rate that the phase of Vy( x ) changes is R \ih\xjjj), which generally does not 

4 We can always make an arbitrary dX conserve the normalization by adding a term that 
is proportional to A. 
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equal to the chemical potential ji y the difference is 


A = Re (ipj | ih | i/jj ) — Hj 

^■)+Re (VhfeK*)-/^ 

= < \ R j | ) + Re < Vh |Qj | ) - hf ] , 


= ( 


n _ 9 

V 2 + C + iA 
2 m 


(5.47a) 

(5.47b) 

(5.47c) 


where 


( 2 ) 

3 


— ( F 


pcs 


|^ 2 |^pc S > 


(5.48) 


The difference A vanishes if we have = QjQ k , but this is generally not true. 


Another problem of Eq. (5.40), that it does not conserve the norm of ^(x), can 
be fixed by absorbing the extra factor into the occupation number o y Assuming that 
Qj 7 ^ Q k for j 7 ^ k, we have 5 

= I = 2 ^ (( I )) = \ Im ^ ^I ^ ’ ( 5 ‘ 49 ) 

and thus we have 




27 

h 


/ Im 


fc=i, 


([^j( x )]Vfc( x )) • 


(5.50) 


The particle number-conserving condition, Ylj=i Qj = 0, follows by noticing that 
£,jk — ikj and I m ([^K x )] 2 ^fc( x )) = — hn ([^(x)] 2 , (/^(x)). Note that the occupation 
number o 3 can be calculated using Eq. (4.47c), 


(2 n) dT (s) 

' i = W) a Sj 


(5.51) 


Taking derivative with respect to the PCS parameter s k , we have 

d Qj = (2 n) d 2 T (2 n) dT dT _ 1 , (2) _ x = ^ 

ds fc T <9s,<9s fc T 2 As- ^ 2n ’ j 

°For the degenerate case, we also need to consider the contributions to q,j from ( i/} k | ipj ) 
and ( i/j k | 4 >j ) for k ^ j. 
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where c jk denotes the covariance of the number operators aj-Oj and a\a k . Using this 
condition, we have 




k =1 

The matrix c is not invertible 


(5.53) 


= = ° =*• c(l,l,...,lf = 0, (5.54) 

k =1 k =1 


but we can invert it in the subspace orthogonal to the vector (1,1 ,..., 1 ) T and thus 
solve for Sj. 


Finally, we come to the point that the orbitals do not remain orthogonal under 
the evolution (5.40). To fix this problem, we notice that any ensemble decomposition 
of the form { Ylk=i V'jkyfOk I 0a ; ), j — 1, 2,..., v } gives the same 1RDM, where U is 
a v x v unitary matrix. For a unitary U generated by —iddt/H, where dj k = d* k j and 
dt the infinitesimal parameter, we have 

dt v 

0'(x) = 00x) - ^ tijkVei&ki*) ■ (5.55) 

in V6j k = i, 

Mi 

By properly choosing the Hermitian matrix d, we can make the orbitals orthonormal 
to each other under the time evolution. Consider the following modified equation for 
the orbitals 


ih^ mod) = 


h 2 

—V 2 + V + R j -d j 
2m 


1 x A 

0i + Qj^*j —— ^-fcV^fc0. 

V 0? fc=i, 

Mi 


/c ? 


(5.56) 


where 


Inn?,- = Im (^ \Qj\ 0* ) = — / Im ^[0*(x)] 2 0^(x)') dx (5.57) 

e i k= 1, J 

Mi 

and 


Re dj = Re (0^ | ih | 0 5 -) — /n,- (5.58a) 

= (0j- l-Rjl 0j ) + Re (0j |Qj 1 0j ) - (5.58b) 

= — f f 2|0i( x )| 2 |0fc( x )| 2 + Re ([0i( x )] 2 0fc( x ))^ x ) - • (5.58c) 

Q i V fc= i 4 / 
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Requiring that the orbitals be orthogonal gives 


o=<*iv;r d) >+<^r d) i^> 


(5.59a) 



(5.59b) 


and this means that the Hermitian matrix d can be determined 



(5.60b) 


(5.60a) 


Note that the degenerate case needs extra care, because both Eqs. (5.49) and (5.60) 
are no longer correct. Instead of making the orbitals orthogonal by small corrections, 
one can just diagonalize the evolved 1RDM in the degenerate subspace (this is similar 
to the nondegenerate perturbation theory). 

Equations (5.50) and (5.56) together determine the dynamics of the PCS ansatz; 
these equations are only valid, however, for cases where v is fixed. This problem 
is caused by replacing the condition (5.30) with (5.31). While the first condition 
can determine | \F err ) in the whole space, the latter one is trivially satisfied in the 
nullspace of the 1RDM and gives no information about | !^ err ) in that subspace. To 
make Eq. (5.31) work in the whole space, we give an infinitesimally small occupation 
number to the states in the nullspace. Thus consider an extra orbital labeled by u+ 1; 
using Eqs. (4.45a) and (4.47c), we have 


<9 2 T 





JO «/|j7|=sin0 


(5.61) 


Supposing that the occupation number of the extra orbital is small, g u+1 <C N, we 
have s u+ i <C Sj for j — 1,2,... ,v. The main contribution to the integral in Eq. (5.61) 
thus comes from 9 ~ 7 t/2 , because of the exponential, i.e., \y\ ~ 1, and we have the 
following approximate relation, 



(5.62) 
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Using this relation and the marginal condition of the 2RDM, we have 


S1/+1 ,j — Pv-t-i Qj 1 


(5.63) 


for j = 1,2,... ,v. The evolution of the occupation number g u+1 thus takes the form 

2 9 


Qu+i 


h 

3= 

^9Qv+i 

h 


:+i(x)] Vb(x) dx 


^+i.i / Im 

.7=1 •’ 

V n 

^jQj / Im ([VC+i( x )]Vj(x)) dx . 

. 7=1 J 


dx 

(5.64a) 

dx . 

(5.64b) 


From Eq. (5.64b) we see that the occupation number g v+l grows exponentially for 
short time, but its initial value must be nonzero to have a nontrivial result. 6 Since 
we have a factorized form (5.63), the chemical potential of the extra orbital is 


hu+i ~ ^ Vv+i 


r 2 

~2in V +V + R ^ 


Vb+i)+ Re (Vb+i |Q^+i| VC+i) , (5.65) 


where 

V V 

i?i,+i(x) ~ 2 Qj l^jWI 2 1 Qu+ i(x) - g ( x ) • ( 5 ‘ 66 ) 

i=i i=i 

A necessary condition for g u+1 to grow steadily is that p u+1 ~ pj for at least one 
j G {1,2,..., u}] otherwise the sign of the integral in Eq. (5.64b) changes rapidly. 


6 The nonzero occupation number can come from two sources: (i) Q v+ \ is small but 
nonvanishing for the PCS ground state; (ii) by going to the next level of approximation of 
Y(s), we can have a nonzero g u+ 1 ~ N°, but N must be finite for the extra level to make a 
difference. 
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Appendix A 

Phase Diffusion: Many GPEs vs 
Bogoliubov 


Good physicists think in many different ways so that they can pick the 
most efficient one. 

- Carlton M. Caves 

The relative phase of the two components of a BEC can be well defined when there 
is an uncertainty in the relative number of bosons. The many-body wavefunction 
thus is a superposition of states with different numbers of particles distributed in 
the two components; such states evolve at different rates due to their different, 
chemical potentials caused by the particle-particle interactions. This effect, called 
phase diffusion, degrades the relative phase of the two components. One approach 
to analyzing phase diffusion is by using many Gross-Pit.aevskii equations (rnany- 
GPEs) [LTRS09], while another approach is to use the number-conserving Bogoliubov 
approximation [Spr02]. Very similar results are derived using these two approaches. 
Here, 1 discuss the equivalence of the two methods using the results in Chap. 3. 

In the many-GPEs approach, we consider the state vector of a two-component, 
condensate of N bosons, with no depletion out, of the two modes, expanded in the 
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Fock basis, 


8>) = V . 



N 

vac) = Xj I j, N - j) , 
i=o 


(A.l) 


where ^ is the creation operator for the crth hyperfine level with the spatial 
wavefunction 0 CT , and the y^s are the complex coefficients that determine the state. 
When the coupling between the two modes in the Hamiltonian (3.69) goes to zero, 
i.e., u 12 = 0, each number sector N l , N 2 ) is decoupled from other sectors. Thus, 
one can solve the whole dynamics by considering the dynamics within each sector 
separately. Because the interaction term is nonlinear in the number of particles, each 
sector evolves according to a different GP equation, 

= (Hi + hw n + gnN^cj)^ 2 + g u N 2 |0 2 | 2 )0i j (A.2) 

ih(f> 2 = i^H 2 + hu 22 + g 22 N 2 10 2 1 2 + ( 721 ^ 1 101102 j (A.3) 


where Hi 2 are the single-particle Hamiltonians. This dependence of GPE on Ni and 
N 2 eventually causes phase diffusion. 


In the Bogoliubov approach, we consider a two-component. BEC of N bosons 
condensed in the wavefunction 

I 0(f)) = ^ (<*i(f) |0i(0)) ® 11) +a 2 (i) |0 2 (*)) ® I 2) j , (A.4a) 

which solves the Gross-Pitaevskii (GP) equations 

+ (A ' 5) 

Here, the GP Hamiltonian takes the form 
(h^ hw 12 \ 

H gP =\ =HS\l)(l\+H$\2)(2\+hw 12 \l)(2\+hw 21 \2)(l\, (A.6) 

\hw 21 Hgp ) 

where 


Hgp — H\ + hw u + 5 , n|tt 1 | 2 |0 1 | 2 + g-[ 2 1 ol 2 | 2 1 02 1 2 , (A.7) 

Hgp = H 2 + hw 22 + ^22 1 <^21 2 1021 2 + 5 , 2ll Q hM0ll” 


(A.8) 
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For the case cu 12 = 0, a 1 and a 2 do not change in time, and the following state also 
solves the GPE (A.5), 

I 0(f)) = ^7 ( “2 I 0i(f)) ® 11) - M I 0 2 (f)) ® I 2 ) ) . (A. 9) 

Recall that the number-conserving Bogoliubov Hamiltonian matrix (3.96) takes the 
form 

(H„ p + \a\ 2 Q$G<h*Q a 2 Q$G<S>Q* \ 

H nch = , (A. 10) 

\ (a) 2 Q*$*G<h*Q H* p + \a\ 2 Q*$*G§Q*J 
where Q(t) = 1 — | 0(f))(0(f) |, and 


<E = - 1 

f cy,\(j)\ 

0 ) 

| . G= 1 

( 9n 

dl2| 

a | 

\ 0 

a 202 ) 

1 1 

921 

922 ) 


(A.ll) 


To separate the part that causes phase diffusion, we divide the Bogoliubov Hamiltonian 
matrix into three pieces, 


H neb ^ncb-lL T ^ ss T Hpd ■ 

The part H nch± is orthogonal to the subspace spanned by | 0) and | 0), 

a 2 mG$R* 


(A.12) 


H 


H gp + \a\ R&G& R 


ncbJL 


(a) 2 R*$*G$*R H* p + |a| 2 /r<h*G<h/r 


(A.13) 


where R(t) = 1 — 10(f))(0(f) | — | 0(f))(0(f) | is the projection operator on the 
orthogonal space. The second term H ss in Eq. (A.12), the spin squeezing term, only 
contains contributions from the mode 0 

|a| 2 | 0)(0 |<EG<E*| 0>(0| a 2 |0)(0|<EG<f>|0*)(0*| 


H aa = 


(ct*) 2 1 0*)(0* |<E*G < 1>*| 0) (01 |a| 2 1 0*)(0* |$*G4>| 0*)(0 : 


2 I / * \ / / * I A 


|a| 2 |0)(0| ae 2 ' l ' e \ 0)(0 


= V 


(a*) 2 e~ 2id | 0*)(0 | M 2 |0*)(0 : 


2 | / * \ / / * 


(A.14a) 


(A. 14b) 


where 6 = arg(aqa 2 /« 2 ), and 

fj = (0 |<EG<E*| 0) = e~ 2ie (0 |$G<h| 0* > = e 2ie (0* |<FG<F| 0) 

I«i| 2 |« 2 | 2 '• 


dll |011 4 + ^22 1021 4 — 2(yf 12 |0 1 | 2 |0 2 | 2 G?X . 


a 


(A.15a) 
(A.15b) 
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The third term H pd in Eq. (A. 12), describing the coupling of the mode </> to the 
orthogonal modes, is 


/ |a| 2 U)(0|$G<r J R a 2 | <j>)(<j> I R* \ 

Hpd = 1 + H.c. , (A.16) 

y(a*) 2 | </>*)(</>* |a| 2 | </>*)(</>* \$*G$R*J 


and the corresponding Hamiltonian reads 


+ (a*)\( 4>* |<T |t() x ) + |a| 2 a 0 ( 0* ]$*G$| ^ ) + H.c. 

= (, ae i6 a\ + ae~ ie a^ (a*e^ e ( 4> |$G«T | t\> ± ) + ae ie (0* |$*G<1>| ) 


(A.17a) 


(A. 17b) 


where i])j_ denotes the field operator orthogonal to the two modes |</>) and \4>). 
Another way to think about the Hamiltonian (A. 17b) is to recall that it arises, in 
the Bogoliubov approximation, from replacing a^ and aj, by a and a* in the original 
Schrodinger-picture Hamiltonian. Restoring the creation and annihilation operators 
for the condensate mode to the Hamiltonian (A. 17) gives 


H p d = (e“’oj,a* + e “’op*) (e “ > op0|4G4*|i))j L ) + e' , <J*(0|$G$|'l4)) . 

(A.18) 

As an example of how the Hamiltonian (A. 18) works, consider the situation where 
the condensate wavefunction </> is an equal superposition of the two hyperfine levels, 


.e., aq = a 2 = 1 /-\/2, 6 = 0, and 


a. 


0 ( 0l H ° 2 ) ’ “'<? ^2 


CLj 


- /., ("I - 02 ) . 


(A. 19) 


where dp (a 2 ) is a shorthand for aq ^ (a 2)< * 2 ). By introducing the Schwinger operator 


Jz = ^ (4 a i - 4 ^ 2 ) = ^ (4^ + aW) > 

we bring the Hamiltonian (A. 18) into the form 


(A.20) 


ftp* = 2 J z (op 0 IKG4* I ih.) + o*< 4 , I*Gt 1 1 ).;)) . 


(A.21) 
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Now we are in a position to see how the many-GPEs approach is related to the 
Bogoliuhov approach. Note that J z counts the difference between the numbers of 
particles in the two components. Thus, we can regard Eq. (A.21) as a controlled 
Hamiltonian depending on the value of N 1 — N 2 , which accounts for the different 
GPEs (A.l) for different sectors. 
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Appendix B 


Bargmann-Fock Representation 
and Its Applications in BECs 


God used beautiful mathematics in creating the world. 

- Paul Dirac 

In Chap. 3, we have introduced the extended catalytic state as a coherent state for the 
condensate mode and a state for the orthogonal modes. The physical state of the BEC 
is encoded in the A-particle sector of the Extended Catalytic State (ECS) Eq. (3.5). 
To decode, one simply projects the extended catalytic state to the IV-particle sector. 
In this appendix, I discuss how to represent the ECS and to project it into number 
states with the Bargmann-Fock (BF) representation [Bar61, Bar62], 

The BF representation, introduced to quantum optics by Glauber [Gla63] , 1 is a 
representation where quantum states of bosonic modes are mapped onto analytic 
(holomorphic) functions, which are convergent with respect to the Gaussian measure 

11 ,/11 gauss : = J l/(-)| 2 e" N d 2 z . (B.l) 

1 The Bargmann representation also finds various applications in fields such as quantum 

foundations [Kem94], quantum gravity [TW01, LT12], quantum chemistry [Har94, MFL06], 
nuclear physics [FL05], and lattice quantum mechanics [CDD + 95]. 
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Many powerful results of analytic functions in the complex plane can thus be exploited 
within the quantum-mechanical context. The BF representation is so versatile because 
we can manipulate it by multiplying functions, taking derivatives, and doing integral 
transformations. Although it is used extensively by mathematical physicists, the 
BF representation has not attracted as much attention from the quantum optics 
and condensed-matter communities, with two possible reasons being: (i) it is not 
very convenient for mixed states; (ii) it is not as intuitive as the quasi-probability 
distributions. We will show, however, that the Bargmann-Fock representation is ideal 
for certain tasks, including many of ours, for which the quasi-probability distributions 
are not well suited. 


In the BF representation representation, the Fock states (number states) are 
mapped to monomials, 


\n i 


,n 2 ,...,n v ) ->■ 1 f 




n k 


and consequently the annihilation and creation operators take the form 


(B.2) 


d 


A 


a k —> -t— , a k —> z k . 

ciZk 


(B.3) 


Since any state in the Fock space can be expanded in the basis of number states, 

I &) = /m,Ba.«„ I ™i, n 2, . , n u ) , (B.4) 

bn-} 


we have the BF representation for that state, 

( v n k \ 

fn x ,n 2 ,...,n v r —f j 

*=1 V 

and the state is reconstructed from the BF representation by 


(B.5) 


&) = <?>(«!, 4, • • • ,al) l vac ) 


t' 


(B.6) 


Another quite useful definition of the BF representation, in terms of coherent 

states, is 


®iip)0i,^ 2, ■■ ■ ,-z») = e |z| /2 . ,z‘„ | V) , 


(B.T) 
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where ( zf, z 2 ,..., z* | is the bra for the coherent state that has complex amplitudes 
(zl, z 2 , ■ ■ ■, z* v ), and |z| 2 = \z 1 \ 2 + \z 2 \ 2 + ■ ■ ■ + \z v \ 2 . Note that the BF representation 
should not be confused with the square root of the Husimi Q function despite their 
similar appearance. Note also that 

e |z| /2 (F\ z k ,z 2 ,...,z*) = [ ( B^ ) (z 1 ,z 2 ,...,z u ]}* = ©**>(*!, z 2 ,..., z v ) , (B.8) 


which defines what we mean by complex conjugation of the BF representation. 


Using the definition (B.7) and the properties of coherent states, it is not hard 
to check Eq. (B.3). Also, using Eq. (B.7) and the completeness condition of the 
coherence states, 2 we have the following result for the inner product of two quantum 
states, 

<!?2 I !? 1 ) = f(V 2 \ z lz' 2 ,...,z;)(zt,z' 2 ,...,z:\VOd\---d\ (B.9a) 

= F t •B|„ a) (z)S| 4 , i ,(z)e _|z| d 2 z . (B.9b) 

Thus the inner product of | 'Tj ) and | F 2 ) is given in terms of the BF representation 
by Eq. (B.9b); indeed, it is easy to verify, by plugging Eq. (B.5) into Eq. (B.9b), that 
the inner product reduces to the standard Fock-space form. 

Consider the operator C = exp ( G* k aja fc ), which is not necessarily Hermitian, 
i.e., G is not necessarily anti-Hermitian. This operator imposes a linear transformation 
on the amplitudes of a coherent state, 


C I 


\ (z T L T L* 

) = e K 


-|z| 2 )/2| 


). 


(B.10) 


G • 

where L = e . Notice that we can obtain any invertible transformation of the 
coherent-state amplitudes by this method. We have 


<^im> 


~u JWCIWW*. 

J *S|^ 2 >(^z) *S|^>(z) e“ |z| d 2 z. 


— j | a) ( a | d“ 


a = — 


a* )( a* | d 2 a = 1 


(B.lla) 

(B.llb) 


2 


7 r 


7r 
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For the same operator £, we also have (z* \C = m 
quently, 


/ l T * T l * 

(z L L z — 


(L 1 z* |, and conse- 



(B.12b) 


(B.12a) 


Thus, we have 


7T 


1 




/ 



(B.13) 


which says that doing a linear transformation on one of the BF variables in the inner 
product (B.9b) is equivalent to doing the conjugate transformation on the other one. 

Many useful properties of the BF representation are discussed in [VB94, Vou06]; 
these include the relation of 93 1 to the wavefunction in the coordinate or momentum 
representation and also the operator-representation in the Bargmann-Fock space and 
its relation to the Glauber-Sudarshan P function, the Husimi Q function, and 
the Wigner W representation. In this dissertation, particular emphasis is put on 
representing Gaussian states and Gaussian unitaries. We also discuss how to project 
an arbitrary pure state to the IV-particle sector with the BF representation. 

Pure Gaussian states have a particularly simple form in the BF representation. 
For example, we have the following for the coherent state | a) = T>(a)\ vac), where 


V(a) = e 


ota, —a a 



—a a 


(B.14) 


is the displacement operator: 



(B.15a) 

(B.15b) 

(B.15c) 


For squeezed vacuum states, 5 ( 7 )| vac), where 


5(7) = e 2<cia 7 “ \ 7 real, 


(B.16) 
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is the single-mode squeeze operator, we have 


®sqv(7>2) = ezN (z*|«S(7)|vac) 


,§H 2 


yj cosh 7 

1 

yj cosh 7 


( I exp (- 


tanh 7 


j 2 j | vac ) 


exp 


tanh 7 2 
-T- Z 


(B.17a) 

(B.17b) 

(B.17c) 


Here we use the following “quasi-normal-ordered” factored form of the squeeze 
operator [Per77, Hol79], 


5(7) = exp (- tai ^ h7 a f2 ) (cosh 7 ) 0 a exp ^*^ 7 a 2 ) . (B.18) 

Notice that the result in Eq. (B.17) can also be derived by Ending a differential 
equation for <B sqv ( 7 ,z): 



270 a vac) = 0 

(B.19a) 

. 1 2 1 f2 

(a cosh 7 + a T sinh 7 ) e 2 “ 27(1 vac) = 0 

(B.19b) 

(^7 cosh 7 + 0 sinh 7 ) <B sqv ( 7 , z) = 0 

(B.19c) 

®sqv(7>^) oc exp ( 2 

(B.19d) 


The magnitude of the proportionality factor can be determined from the normalization 
condition 


IJ l*£ sq v(7,-)| 2 e N d 2 z= 1, (B.20) 

and the phase of the proportionality factor can be fixed by knowing that *B sqv ( 7 , 0 ) = 
l/y'coshy is real and positive. 


It turns out that the BF representation of the most general pure Gaussian states 
for an arbitrary number of modes (not at all restricted just to coherent and squeezed 
vacuum states) can be determined from the action of Gaussian unitaries on the 
vacuum. Gaussian unitaries are generated by Hamiltonians containing terms linear 
or quadratic in the annihilation and creation operators; they take simple forms in the 
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BF representation. For the displacement operation (B.14), we have 
displaced by a 




®2}(a)| ijj) \Z) 

— i|a| 2 olz —a*d/dz <vv / \ 

= e 21 1 e e ' *B| ^ ) {z) 

®coh(^b %) Oi ) ■ 


The rotation (phase-shift) operation, 


77(0) = e 


id a 1 a 


is also straightforward: 

, . rotated by 0 _ ie 

^li>)K z ) ^ ®72(0)| ip) GJ ®| ip)\^ %) • 


(B.21a) 

(B.21b) 

(B.21c) 


(B.22) 


(B.23) 


This result can be derived easily by using the definition (B.7). For the squeezing 
operation (B.16), we have 

®|„M ®« t)W M (B.24a) 

= 6XP (i £2 -\ z ‘ 1 ')^)( z ) (B.24b) 

= ®sqv(7,-) (coshy) - ^ exp^ ta1 ^ 17 ®^)(z) , 

(B.24c) 

where the factored form (B.18) of the squeeze operator is used. The term (cosh y) _zd//dz 
rescales the independent variable as z —> z/ coshy; the general rule here is that for r 
real, r a a is represented by T zd / d ~ i n the BF representation and that 

T zd/d ‘f(z) = }(zt) . (B.25) 

The other term in Eq. (B.24), exp(^ tanhq d 2 /dzj , can be evaluated in the Fourier 
domain; it maps Gaussian functions to Gaussian functions without changing their 
centers. 


Another important Gaussian unitary operation is the passive linear-optical network 

or multiport beamsplitter. The beamsplitter is described by a unitary operator LA, 

which acts on v input modes according to 

multiport beamsplitter + v - ' rr rn r, r \ 

a k -----» Wa k U — } y cij U jk . (B.26) 
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Here U is a matrix that must be unitary in order to preserve the canonical commutation 
relations. This transformation can be written in the equivalent, inverted form 

V 

Ua{U ] = Y,U kj a). (B.27) 

3 = 1 

Using Eq. (B.6) or Eq. (B.7), we can find the transformation law for the BF represen¬ 
tation, 


-,Z V ) 


multiport beamsplitter , . 

= ®| z' 2 , ■ ■ ■, z' u ) , 


(B.28a) 

(B.28b) 


where 


'y ^ U kj Zj . 

3 


(B.29) 


Most generally, Gaussian unitaries are unitaries generated by Hamiltonians con¬ 
taining terms linear or quadratic in the annihilation and creation operators. It turns 
out that the simple transformation rules of displacement in Eq. (B.21), single-mode 
squeezing in Eq. (B.24), and beamsplitters in Eq. (B.28) are enough to construct 
the most general Gaussian unitaries. Using displacement operations, we can remove 
the linear terms in the generators of a Gaussian unitary, and the resulting unitary is 
generated by a quadratic Hamiltonian. The Bloch-Messiah reduction theorem [BM62] 
states that any multimode Gaussian unitary generated by a quadratic (Bogoliubov) 
Hamiltonian can always be decomposed into a multiport beamsplitter, followed by 
a set of single-mode squeezers on each mode, followed by yet another multiport 
beamsplitter, 

W bog = u msp f n S k('rS) VLp • (B.30) 

' k =1 ' 

The BF representation allows a simple manipulation of all these elementary operations 
and thus is a suitable platform to work on Gaussian transformations. 

Any pure Gaussian state can be generated by a Gaussian unitary acting on the 
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Multiport Single-mode Multiport 

splitter squeezers splitter 


Mode 1 


Mode 2 


Mode u 



Figure B.l: The Bloch-Messiah reduction theorem states that any Gaussian 
unitary generated by a quadratic Hamiltonian can be decomposed into a multiport 
beamsplitter V , followed by a set of single-mode squeezers on each mode, followed 
by yet another multiport beamsplitter U. 


vacuum state," 


| F) =U | vac) (B.31a) 

= £>(a 1; a 2 , • • •, ay) U hog | vac) , (B.31b) 

where the Gaussian unitary is broken up into a displacement operator and a unitary 
lA hog which is generated by a quadratic (Bogoliubov) Hamiltonian. Using the Bloch- 
Messiah reduction theorem, we have 


F ) = V(a 1 ,a 2 , ...,oi v )U, 


msp 


0^2? • • • 5 


msp 


n s k { 7 fc 

k =1 
v 

n 


V 


k =1 


'Lp I vac) 
vac ) . 


(B.32a) 

(B.32b) 


where the initial multiport beamsplitter Vn lsp can be discarded because it does 
not change the vacuum state. All the elementary operations in Eq. (B.32b) have 


^The covariance matrix of any pure Gaussian state is symplectically equivalent to that of 
the vacuum state, and this symplectic transformation can always be realized by a Gaussian 
unitary. 
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been discussed before, in Eqs. (B.17), (B.21), and (B.28). Thus we have the BF 
representation for the most general Gaussian state, 


&\9)(Zi,Z 2,-.-,Z u ) = 


®coh(^fe; %k) 

L J[ V cosh 7^ 


exp 


tanh 7 fc 


Ukj ( z j 01 j) 

3 =i 


(B.33) 


For the Extended Catalytic State (ECS) Eq. (3.5), where the state of the noncon¬ 
densate modes is a squeezed state centered at the origin, we have 

®ecs(^0> ■ 1 “u) 


®coh(®> ^o) 


L = i Vcoshqfc 


exp 


tanh r ) k 


Ukj z ,j 

3 = 1 


(B.34) 


where the condensate mode is denoted by j = 0 (and hence with the BF variable 
z Q ) and the other modes, orthogonal to the condensate mode, by j = 1, 2,..., u. 
Projecting the state | F ecs ) to the IV-particle sector is equivalent to keeping only the 


terms of power N in the BF variables in a power-series expansion of *B ecs (z). 


With the BF representation, it is also convenient to calculate the correlation 
matrices or, equivalently, the reduced density matrices (RDMs) of a multi-boson 
system. For example, the single-particle reduced matrix (1RDM) takes the form 


Pkj = ('P\ a ] a k\'P) (B.35a) 

~ 3jk T J ®jtf>(z)9V)(z )zlz j e~ lzl d 2 z. (B.35b) 

Returning to the case of the extended catalytic state, where the physical state is 
encoded in the IV-particle sector of the extended catalytic state, we have 

Pkj = (^ecs \V N a)a k V N | R ecs ) (B.36a) 

= (^ecs I a)a k V N \ R ecs ) . (B.36b) 

A convenient way to evaluate Eq. (B.36b) is through the generating function 

G k ,j(r) = ( ^ ecs | a)a k r M \ F ecs ) (B.37a) 

= (^ecs | t M/ 2 a)a k r M/2 1 F ecs ) (B.37b) 

= J l^ecslzv^r) | 2 (zjZ* k - 5 jk ) e _|z| d 2 z , (B.37c) 
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where A f = JT a ] a j I s the total particle-number operator, so that r A = Y1n=i t n V n . 
The value of [R- is encoded in the coefficient of the term t n in a power-series 
expansion of G k j(r). For Gaussian states, the integral (B.37c) is not hard to do, 
thus allowing to be retrieved. 

This procedure of deriving lRDMs by using the BF representation can be straight¬ 
forwardly generalized to the higher order RDMs. 
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Appendix C 

Symplectic Methods for Quantum 
Mechanics 


Mathematicians must have felt this way when they discovered that com¬ 
plex numbers were more than just one extra gimmick: Virtually every 
idea of mathematics, from the geometry of curves to the analysis of par¬ 
tial differential equations, was ripe for complexihcation. Mathematics 
exploded overnight. 

- fan Stewart 


C.l Introduction 

Quantum held theory is notoriously hard, which explains why superconductivity 
and superfluidity eluded people for so long. Our knowledge of them, however, has 
tremendously increased with the help of the Bogoliubov transformation [Bog47, 
BCS57b, Kit87] . The Bogoliubov transformation of n bosonic modes is a linear 
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transformation of the annihilation and creation operators, 1 

a = Ub + V*b\ (C.la) 

= U*b ] + Vb , (C.lb) 

where a = (% a 2 • ■ ■ a n ) T and — (a\ a\ ■ ■ ■ a} n ) T are the column vectors of the 
input annihilation and creation operators, b = (bi b 2 ■ ■ ■ b n ) T and = (b\b\ ■ ■ ■ bh) T 
are the output operators determined by implicitly Eqs. (C.l), and U and V are n x n 
matrices that specify the transformation. The inputs, being annihilation and creation 


operators of different modes, satisfy the canonical commutation relations, 

[®j, {P'j i 0) \^j i ^jk • (C.2) 

The Bogoliubov transformation Eqs. (C.l) conserves these relations, 

[6j,6*] = [6;,6i]=0, [&„ &1] = ^ , (C.3) 

or equivalently, we have the following constrains on the matrices U and V, 

U f U-V f V = 1, (C.4a) 

U T V - V T U = 0 , (C.4b) 


where 1 and 0 denote the n x n identity matrix and null matrix, respectively. These 
conditions (C.4a, C.4b) are called the symplectic conditions. Therefore, the study of 
Bogoliubov transformations is equivalent to the study of symplectic transformations. 

Symplectic structures arise naturally in classical Hamiltonian mechanics, where 
linear transformations of the canonical coordinates conserve the Poisson brackets. Not 
restricted to Hamiltonian mechanics, it also finds applications in many other aspects 
of physics and engineering; the reader is referred to [BerOO, Arn89] for pedagogical 
introductions to the subject, [FOP05, Gos06, Xia09, ARL14] for applications to 
quantum mechanics, and [GS90] for general applications to physics. 

1 For the purpose of this dissertation, I only discuss the bosonic case and leave out the 
equally important fermionic case. 
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This Appendix provides a systematic introduction to the symplectic formalism to 
quantum mechanics in the basis of the annihilation and creation operators. 2 Doing 
this is especially convenient for the purpose of quantum optics and quantum rnany- 
body physics, partly because the many-body Hamiltonians are usually written in 
terms of the annihilation and creation operators. Another, perhaps more important 
reason is that using the basis of annihilation and creation operators allows one to 
distinguish beamsplitters from squeezers easily; while beamsplitters are considered 
readily available in quantum optics, squeezers are much harder to implement and 
considered as a resource [Bra05]. 

Following is a list of the topics covered in this Appendix: (i) the standard symplec¬ 
tic structure of a quadratic Hamiltonian; (ii) the Heisenberg time evolution under a 
quadratic Hamiltonian; (iii) the relation between the “diagonalization” of a quadratic 
Hamiltonian and the symplectic eigenvalue problem; (iv) the polar decomposition of 
symplectic matrices with a proof of the Bloch-Messiah theorem [BM62]; and (v) some 
examples and applications to quantum optics and quantum many-body physics. 


C.2 Quadratic Hamiltonians 


The class of Hamiltonians that are quadratic in annihilation and creation operators, 
or simply quadratic Hamiltonian, is a very special one. Many important Hamiltonians 
can be reduced to this class if approximations are allowed. Most importantly, systems 
evolving under quadratic Hamiltonians can be solved efficiently, and their ground 


states are easy to find. Before going on, however, let us look at the notation for an 
n-rnode quadratic Hamiltonian, 3 


n 





(C.5) 


2 

Almost all materials on symplectic structures adopt the position-momentum basis of x 
and p. 

3 

For quadratic Hamiltonians, the difference between different orderings of the annihilation 
and creation operators is a c-number. In this Appendix, we assume all Hamiltonians are 
symmetrically ordered. 
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where H x , H 2 , H 3 , and are all n x n matrices, and the following conventions are 
used throughout this Appendix, 



\al 

The 2 n x 2 n matrix in Eq. (C.5), often called the Hamiltonian matrix, determines 
the quadratic Hamiltonian TL. To distinguish it from ordinary n x n matrices, we 
denote it by a special sans serif font, 4 


H 




(C.7) 


The canonical commutation relations (C.2) take the form 



where the commutator notation on the left means to put the commutators of the 
elements in the column and row vectors into the corresponding positions of the 2 n x 2 n 
matrix. We also have 



\ ®ra ) 


(C.9a) 


(C.9b) 


4 In this dissertation, slanted sans serif fonts always denote matrices of symplectic 
structure. 
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where 

/ 0 1 \ 

X= . (C.10) 

0 j 

In words, Hermitian conjugation of the annihilation and creation operators is the 
same as multiplication by X. Equations (C.8) and (C.9) introduce the matrices 
Z = = Z* = Z 1 and matrices X = X^ = X* = X 1 in the ways that demonstrate 

their importance within the symplectic formalism. 

Using Eqs. (C.9), it is easy to see that matrix XH 1 X gives the same Hamiltonian 
hi as does H. Thus we can redefine the Hamiltonian matrix by 

H^^(H + XH t X). (C.ll) 

The new H, which gives the same Hamiltonian hi. is now constrained by the condition 


H t = XHX. 


(C.12) 


With this condition, the Hamiltonian matrix H has a one-to-one correspondence with 
the quadratic Hamiltonian hi. Given the condition (C.12), the Hermiticity of hi 
requires H to be a Hermitian matrix, 

H ] = H. (C.13) 

Combining Eqs. (C.12) and (C.13), we have 

H* = XHX. (C.14) 


Only two of the three conditions, (C.12), (C.13), and (C.14), are independent; we 
use the conditions (C.13) and (C.14) as primary (the reason for this choice should 
become apparent after the introduction of symplectic matrices). Imposing these two 
consraints, the Hamiltonian matrix H adopts the form 


H = 


h x m 


h 2 hi 


where H 1 = // j. and H 2 = H 2 . 


(C.15) 
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The Hamiltonian hi therefore takes the form 


« = S 





(C.16) 


h x h; 
h 2 m t 

It is clear that H 1 describes an n-port beamsplitter and H 2 a squeezer. In contrast, 
had we written the Hamiltonian in the basis of x and p, it would not be so easy 
to distinguish these two different sorts of transformation. We have now prepared 
ourselves with the representation of the quadratic Hamiltonian hi, and the next step 
is to solve for the time evolution. 


C.3 Heisenberg Time Evolution 


In the Heisenberg picture, the annihilation and creation operators undergo a linear 
transformation under the evolution of the quadratic Hamiltonian hi(t), 




(C.17) 


where S(t) is the evolution matrix. Note that the initial values are used if the time 
arguments of the annihilation and creation operators are omitted. Because any 
operator on the Fock space can be expanded as a Taylor series of the annihilation and 
creation operators, knowing S(t) allows one to solve the whole dynamical problem.' 


Consider the following time-dependent quadratic Hamiltonian in the Schrodinger 
picture 



(C.18) 


Going to the Heisenberg picture, we have 
1 


^o(t) = g ( aT W a C0) 


H,(t) H*{t) 
H 2 (t) m(t) 



(C.19) 


"’Different ordering of the annihilation and creation operators results in different expansion 
coefficients. 
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Since the annihilation and creation operators in the Heisenberg picture also satisfy 
the canonical commutation relations, we have 


[n o {t), a(t)] = —H^t) a(t) - H* 2 (t)a\t), (C.20a) 

[H 0 (t), a f (t)] = H 2 (t) a (t) + m(t) a \t) . (C.20b) 


Thus, the Heisenberg evolutions of the annihilation and creation operators obey the 
differential equations (H = 1) 


V at W/ 

Comparing Eqs. (C.17) and (C.21), we have 


(C.21) 


S(t) = - iZH(t)S(t ) , 

which can be solved formally if H is time independent, 6 


(C.22) 


5(t) = exp (-itZH) . 


(C.23) 


The evolution matrix has a symmetry that comes from the fact that a (t) and 
a'(t) are Hermitian conjugates of one another: 


S*(t) = XS(t)X. 


(C.24) 


More intuitively, this symmetry can be represented by the following matrix form, 


S{t) 


fs^t) s 2 *(f)\ 

\S 2 (t) 


(C.25) 


where S 1 (t) and S 2 (t) are n x n matrices. The careful reader will hnd that this 
condition is similar to condition (C.14) on /-/, and it can actually be derived from 
that condition. In the time-independent case, the derivation is 


S*(f) = exp ( itZH*) = exp {itZXHX) = exp (- itXZHX ) = XS(t)X . (C.26) 


6 The formal solution involves a time-ordered exponential when the Hamiltonian matrix 
H(t) is explicitly time dependent. 
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This explains why we choose Eq. (C.14) over Eq. (C.12): the evolution matrix S[t) 
does not satisfy a condition similar to Eq. (C.12). 


In the Heisenberg picture, the canonical commutation relations take the form 


a (t) 
a \t), 


, (a f (f) a(t)) 


= 5(t) 




a^ a 


S\t ) = S(t)ZS\t ) . (C.27) 


Conservation of the canonical commutations means that 


5(t)ZS\t) = Z . 


(C.28) 


This, called the symplectic condition, is a consequence of the Hermiticity of H, 


S(t)ZS\t) = exp (—itZH) Z exp (itH^ Z^) (C.29a) 

= exp (— itZH) exp ( itZH ) Z = Z . (C.29b) 

Two different, but also very useful, ways of writing the symplectic condition are 

5~\t) = ZS t (f)Z, (C.30) 

S\t)ZS(t) = Z. (C.31) 


Note that the determinant of a symplectic matrix is either 1 or —1 as a consequence 
of Eqs. (C.24) and (C.28). 

We will see that the Bogoliubov transformation is a symplectic transformation of 
the modal operators and can be generated in the same way as 5{t) in Eq. (C.23). 


C.4 Bogoliubov Transformations 

We have already seen how to put a quadratic Hamiltonian into a standard form and 
how formally to derive the Heisenberg dynamics of the annihilation and creation 
operators. Questions still remain: (i) how to solve the energy spectrum of the 
Hamiltonian hi and (ii) how to derive S(t) beyond the formal result in Eq. (C.23). 
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To answer these questions, we will need to “diagonalize” the Hamiltonian matrix H 
by the Bogoliubov transformation (C.l), 


= B 



where B = 


U V* 
V U* 


(C.32) 


The matrix B satisfies the condition 


B* = XBX , (C.33) 

by virtue of the top and bottom halves of the vectors being Hermitian conjugates on 
one another, and the preservation of commutators implies that [this is equivalent the 
condition (C.4)] 

B ] ZB = Z . (C.34) 

Thus the matrix B shares the same properties as S(t) and is also symplectic. 

In the new basis of b and b^, the Hamiltonian matrix becomes 

H' = B ] HB. (C.35) 

The condition (C.14) is preserved by this symplectic transformation, 

XH'X = XB ] HBX = B t XHXB* = H'* . (C.36) 

It is also straightforward to see that H' = H ''. Therefore, symplectic transformations 
preserve the symmetries of the Hamiltonian matrix. 

Often it is useful to consider conserved quantities under symplectic transformations. 
For example, symplectic transformations keeps the determinant of the Hamiltonian 
matrix unchanged, 

det(H') = det (H). (C.37) 


Also, we have 

tr [ZH'ZH') = tr(ZB ] HBZB ] HB) = tr (ZHZH) 


(C.38) 
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which can be easily generalized to higher orders. Note that the trace of any odd 
repetition of ZH is zero, for example tr (ZH) = tr(Z HZHZH) = 0.' We will use these 
invariants in Secs. C.9.1, C.9.4, and C.9.6 to solve practical problems. 

The condition Eq. (C.34) may not be the most convenient for the purpose of 
expressing the symplectic restrictions on £>; instead, one can use the exponential form 

B = e iZM , (C.39) 

where the matrix M has the same constraints as the Hamiltonian matrix (indeed, this 
the same trick as writing a unitary operator U as U = e ~ lH , where H is Hermitian), 

M* = XMX and = M. (C.40) 

We will use this form to parameterize all the two-dimensional symplectic matrices in 
Subsec. C.9.1. 


C.5 Symplectic Form and Basis 

The symplectic form, which plays a role similar to the inner product, is defined as 


™b P {Qi, &) = {Q2\Z\qi), (C.41) 

where | q 1 ) and | q 2 ) are 2n-dimensional vectors in the symplectic vector space. A 
symplectic transformation is a transformation that preserves the symplectic form, 

(92 1^ I } — ( 92 |S f ZS | 9,) — (q 2 \Z \ q 1 ). (C.42) 

Another useful concept is that of a symplectic basis, 

{ I ) | 3 = 1, 2 , ■ ■ ■, n; ° = 0,1} , (C.43) 

which satisfies 

(4 fc) |z|sg ) ) = (-i)x«^,fc- ( c - 44 ) 


( One can simply use H = XHX, H' = H, and XZX 


Z to prove this. 
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For a particular j, the two values of a label two symplectically conjugate states, 
which are related by 

(C.45) 

(cr + 1 should be interpreted as modulo 2 addition; i.e., a changes to the other value). 
An example of a symplectic basis is the natural basis 

\e { J ) ) = (0 0 ■■■ 0 10.0 0) T , (C.46) 

where there is just one nonzero element, equal to 1, located at the (an + j) th position. 
A symplectic transformation is a transformation that transforms a symplectic basis 
to another symplectic basis, 

I ) = S | s ( J ] ). (C.47) 

It is straightforward to verify that Eqs. (C.44) and (C.45) hold for the new basis. 
Note that conjugate states are mapped to conjugate states, 

s | 4i, > = S X I s<?> )* = X (S sf ))* . (C.48) 

Knowing the two bases in Eq. (C.47) allows one to write the symplectic transformation 
matrix as 

n 

S = E (- 1 )” E I X I Z . (C.49) 

cr =0) 1 j =i 

which will be used to find the canonical form of the Hamiltonian matrix in Sec. C.7. 


C.6 Symplectic Eigenvalue Problem 

Knowing the eigenvalues and eigenstates of ZH allows one to solve explicitly for the 
time evolution 

5{t) = exp (— itZH) . (C.50) 

Moreover, these eigenvalues and eigenstates can be used to “diagonalize” the Hamil¬ 
tonian hi. 
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The symplectic eigenvalue equation reads 

ZH \ q) = \ q \q) , (C.51) 

where | q) is an eigenvector and X q its eigenvalue. An equivalent, perhaps more useful 
way to write the eigenvalue equation is 

H\q) = \Z\q). (C.52) 

To find all the eigenvalues, one can solve the characteristic equation 

det (H - AZ) = 0 . (C.53) 

Since H and Z are both Hermitian, we have 

[det(H - AZ)]* = det {H* - A *Z*) = det (H T - A *Z T ) = det (H - A*Z) , (C.54) 

which implies that if A is an eigenvalue, then A* is also an eigenvalue. Using the 
condition Eq. (C.14), we have 

HX\q)* =XH*\q)* = X* XZ | q )* = -X* q ZX\q)* , (C.55) 

which means that —A* is also an eigenvalue of Z H. corresponding to the conjugate 
eigenvector X \ q )*. Thus we have 

H | q) — —A*Z | q) , where | q) — X \ q)* is the conjugate of | q) . (C.56) 

Thus, the eigenvalues appear in quadruplets 

{A, A*, -A*, -A} . (C.57) 

This is an expression of the fact that the characteristic function (C.53) is even and 
real. 

For two eigenvectors | q 1 ) and | q 2 ), we have 

{(h\H\qi) = (q 2 \Z\qi) = A* 2 ( q 2 \Z | ) , 


(C.58) 
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which says that the symplectic form is zero when the eigenvalues of two eigenvectors 
are not a pair of complex conjugate numbers, 

(9 2 |Z | ) = 0 ifA^A*. (C.59) 

A special case is when | q 1 } and | q 2 ) are the same state, and we have 

(q \ Z | q) — 0 if X q is not real. (C.60) 

It turns out that the eigenvalues of Eq. (C.52) can be real, pure imaginary, or 
complex; I will be treat these possibilities separately in the following subsections. The 
central topic is to construct a symplectic basis from the eigenvectors. Such symplectic 
bases turn out to be crucial to bringing the Hamiltonian matrix H to the so-called 
canonical form. 

C.6.1 Real-Eigenvalue Case 

When A is real, the four eigenvalues in the quadruplet Eq. (C.57) degenerate into 
{A, —A}, and the eigenvalue equations take the form 

H \r a ) = (—1) CT A Z \ r a ) , (C.61) 

where a = 0,1 denotes the two conjugate eigenvectors, which are related by the 
condition (C.56), 

I r a+l ) = X \ r a )* . (C.62) 

As a consequence, we have 

( r i|Z|r 0 ) = 0, (C.63a) 

(r 1 \Z\r 1 ) = - (r 0 \Z\r 0 ) . (C.63b) 

Note that the symplectic form of | r a ) with itself, (r a \Z\r a ), must be nonzero 
because of the symplectic nondegenerate argument. 8 Therefore, we can “normalize” 

8 In other words, the eigenvectors defined by Eq.(C.52) span the whole symplectic space. 
Only the null vector has zero symplectic form with all vectors in a symplectic space. 
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the two eigenvectors, so that the symplectic forms in Eqs. (C.63) become standard, 

(v|Z|r ff ) = (-ir<5 ff y (C.64) 

Note that we adopt the convention that the eigenvector that has positive symplectic 
form with itself is labeled as a = 0. The two eigenvectors, | r 0 ) and | r 1 ), form a 
symplectic basis on a two-dimensional plane. 

C.6.2 Pure-Imaginary-Eigenvalue Case 

When A is pure imaginary, the four eigenvalues also degenerate into {A, —A}, and 
the eigenvalue equations are 

H\g ± )=±XZ\g ± ). (C.65) 

Since A = —A*, the eigenvectors are now conjugate to themselves, 9 

\g ± ) = X\g ± )\ (C.66) 

which leads to the following relation 

{(g- \Z\g + ))* = (g_\XZ*X\g + ) = -(g_\Z\g + ). (C.67) 

Thus (g_ \Z\ g + ) is pure imaginary and can always be written in the following 
standard form by rescaling either | g + ) or \g_ ), 10 

(g_\Z\g + ) = i. (C.68) 

Furthermore, from the “orthogonality” condition Eq. (C.60), we have 

(g + \Z\g + ) = (g_\Z\g_) = 0. (C.69) 

9 Strictly speaking, the two sides of Eq. (C.66) can differ by a phase factor e , but it 
always can be eliminated by redefining the eigenvector as | g) —> e~ ld ^ 2 \ g). 

10 Again, we can argue that this symplectic form must be nonzero by the symplectic 
nondegenerate condition. 
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By introducing the new basis 


f 12°) = (1 9 + ) + * 1 9 ~ )) ’ 

(0.70) 

1 ^) = 7/f ( \ g + > - *\ g - >) ’ 


we have 


\9 0 ) = X\9i)*, 

(C.71a) 

(9<r' \Z\ g a ) = (-1 

(C.71b) 

where cr = 0,1. Thus \ g 0 ) and (p ) form a symplectic basis; 

in this basis, the 


eigenvalue equations (C.65) become 


H\g <7 ) — XZ \ g a+1 ). (C.72) 


C.6.3 Complex-Eigenvalue Case 

When A is complex, the eigenvalues appear in a group of four {A, A*, —A*, —A}, and 
the eigenvalue equations read 

H | c a ± ) = ^(—1) CT Re A ± i Im A j Z | c CTj± ) , (C.73) 

where ± denote the complex conjugate eigenvalues and cr = 0,1 labels the conjugate 
eigenvectors. The only situation such that the symplectic form of two eigenvectors 
can be nonzero is when their eigenvalues are a pair of complex conjugate numbers, 

( c„,_ \Z I c 0 , + ) = { d,_ \‘XZX I c li+ >• = - (< d,_ \Z I d, + ) )• . (C.74) 

By multiplying the eigenvectors by appropriate complex numbers, we can always 
“normalize” these symplectic forms, 

(c 0i _ |Z | c 0i+ ) = - (c x _ |Z | c lj+ ) = 1 . 

We introduce the vectors | c a r ) with r = 0,1 as 

^<7,0 ) ( I ^<7,+ ) “h | C(J,— )) / V^2 , 

c<T, i} = (I C CT + 1 ,+ ) - | c ff+li _ ))/V2 


(C.75) 


(C.76) 
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which form a symplectic basis, 

\c a+ l,r) = X\c«, T T , 

( c a >j |Z | c ffiT ) = (-1 ) a 8 a y 5 r y . 

In this basis, the eigenvalue equations (C.73) become 

H | c ffiT ) = (-1) CT+T Re(A) Z | c CT)T ) + i Im(A) Z | c ff+l!T+1 ) . 


(C.77a) 

(C.77b) 


(C.78) 


C.7 Canonical (Normal) Forms 


We have already seen how to construct symplectic bases from the eigenvectors of ZH 
with real, pure imaginary, and complex eigenvalues. These cases can coexist, and the 
whole symplectic vector space is a direct sum of the three. From the “orthogonality” 
condition (C.59), we have the following completeness condition 

( n r n g n c 

E14 >< 4 i+El 4> x4> I+E E 41 ) (41 

j =1 k =1 1=1 r=0,l 

where / denotes the 2 n x 2 n identity matrix, and n r , n g , and n c are the numbers of 
sets of real, pure imaginary, and complex eigenvalues, respectively. Note that these 
numbers must satisfy 


)z, (C.79) 


n r + n g + 2 n c = n . 


(C.80) 


With the completeness condition Eq. (C.79), we have the symplectic spectral decom¬ 
position for the Hamiltonian matrix 


H = H E (-!)" ( E I4’ > < ■I + E19?’ > < 9?> I + E E I ) < 

<7—0,1 ' j= 1 k=1 1=1 r= 0,1 


(,) ''44 I Z 


zee E I E xE I + Ef-DW 1 1 47 xE'i + E E 

<7—0,1 ^ j = 1 fc=l Z=1 T— 0,1 


(-1) T Re A? | c ( JJ T ) (<£ | + {-l) a i Im A ( c ° | c^ 1>T+1 ) (c| 


.(0 


z. 


(C.81) 
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where Eqs. (C.61), (C.72), and (C.78) are used for actions of H on the symplectic basis. 
If ZH has a nilpotent subspace, its eigenvectors may not span the whole symplectic 
space. In this case H cannot be written in the form (C.81) in that subspace; we will 
come to this case at the end of Sec. C.9.1. 


To take the Hamiltonian matrix to the canonical form, we introduce the symplectic 
transformation 


e = E EiTxeg’i+^isfxey^i + Eihlx 


(l) \/_( n r+ n a+0 


<7=0,1 V i=l 


k =1 


1=1 


(C.82) 


_j_ I JO w (n r +n g +n c +l) 

"T" | c <7,l /\ e <7 


where { | ) | j — 1, 2,..., n; a — 0,1 j is the natural basis. The symplectic matrix 

B transforms the basis in Eq. (C.81) into the natural basis, and we have 




(C.83) 


In matrix form, the symplectic matrix B reads 


(u V* 
yc u* 


(C.84) 


where the symmetry XBX = B* is a consequence of the symmetry (C.45) of the 
symplectic basis. 


It is worthwhile writing the canonical Hamiltonian of all the three cases in explicit 
matrix form. For the real-eigenvalue case, 


H = 

''real 



5 


for the pure-imaginary-eigenvalue case, 


= 


'imag 






(C.85a) 


(C.85b) 
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and for the complex-eigenvalue case, we have 


H, 


comp 


/ Re A 
0 
0 


0 0 
— Re A % Im A 
—i Im A Re A 


\—i Im AO 0 

Note that by rephasing the basis set, one 
Eq. (C.85b) real, 


i Im A ^ 
0 
0 

— Re A J 


(C.85c) 


can always make the pure-imaginary form 


H, 


imag 


( 0 |A|\ 

\|A| Oj 


(C.86) 


The only one of the canonical forms (C.85) that can be positive definite is the 
real-eigenvalue case. Because the symplectic transformation Eq. (C.83) keeps the 
positiveness of H, we have 

H > 0 ==> all the eigenvalues of ZH are real. (C.87) 

Thus, any H > 0 can be turned into a collection of independent harmonic oscillators 
with positive frequencies by a Bogoliubov transformation. 


While the real-eigenvalue case corresponds to harmonic oscillators, the pure- 
imaginary-eigenvalue case corresponds to single-mode squeezers. Perhaps the most 
interesting case is when the eigenvalues are complex Eq. (C.85c), 

Hcomp = Re(A) (a* a — b^b) + i Im(A) ( a)tf — ab ) . (C.88) 

In this situation, the Hamiltonian matrix cannot be decoupled into a sum of single¬ 
mode Hamiltonians. Physically, the Hamiltonian (C.88) arises when two sidebands 
symmetrically placed about a carrier frequency, with the carrier energy removed, 
interact via a two-mode squeezing interaction. 


C.8 Polar Decompositions 


The polar decomposition allows one to break a complicated symplectic transformation 
into several simple ones. It can also be used to determine the squeezing ability of a 
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unitary Gaussian operator. Both of these have important applications in quantum 
information processing for continuous variables. 

The polar decomposition of a symplectic matrix S is as follows, 

S = KU = UK (C.89) 

where U and U are unitary, and K and K are Hermitian. Interestingly, the matrices 
U, U, K, and K are also symplectic. The proofs of this for the two ways of polar 
decomposition are similar, so I just sketch the first one here. Since S _1 = ZS^Z, we 
have h/ _1 K" _1 = ZU^ZZK^Z. Because of the uniqueness of the polar decomposition, 
we conclude that U -1 = ZU^Z and = ZK^Z, and these are nothing but the 
symplectic conditions for U and K, respectively. 


The polar decomposition can be used to break down any Gaussian unitary into a 
sequence of multiport beamsplitters and single-mode squeezers. To see this, we first 
consider the unitary part, 

U = exp (- iZH v ) . (C.90) 

Because U is both unitary and symplectic, its generator satisfies 

(ZH.y — ZH V , Hl = H ut Ht = XH v X. (C.91) 

These three conditions imply that H i; must have the following block-diagonalized 
form, 



(C.92) 


where H v is an n x n Hermitian matrix. Thus U is also block diagonalized and 
describes a multiport beamsplitter, 



(C.93) 


For the Hermitian part of the polar decomposition, we have 
K = exp (— iZHff) . 


(C.94) 
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Because K is both Hermitian and symplectic, its generator satisfies 

(■ ZH K y = -ZH K , Hi = H k , H*=XH k X. (C.95) 


These conditions imply that H K must take the form 


H k = 




(C.96) 


rjn 

where H K = H K is an n x n symmetric matrix. Symmetric matrices like H K can 
always be diagonalized by unitary matrices in the following way (see App. E), 

H k = VD V T , (C.97) 


where = V is a unitary matrix and D a diagonal matrix. Thus we have 


H k = VDV' 


where V = 



and D = 



(C.98) 


Putting Eq. (C.98) into Eqs. (C.94) and (C.89), we have 

K = Vexp{-iZD)V ] . (C.99) 

Defining 1/lA = V f U gives us 

S = 1/exp (- iZD ) V ] U=V exp (- iZD ) tV f , (C.100) 


where 1/1/ and V are both unitary and symplectic, and D is of the form in Eq. (C.98). 
The decomposition (C.100) is called the Bloch-Messiah reduction [BM62, Bra05]. 
ft states that a multi-mode Bogoliubov transformation (Gaussian unitary) can be 
decomposed into a multiport beamsplitter, followed by a set of single-mode squeezers 
on each mode, followed by yet another multiport beamsplitter. 

Consider the time evolution (C.22) for the symplectic matrix S, 

S(t) = - iZH(t)S(t ) (C.101) 


where /7(f) is a time-dependent Hamiltonian matrix. Noticing that 
K 2 (t) = S(f)S f (f) 


(C.102) 
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we have 


dK 2 

dt 


55 ] + SS f = - iZHK 2 + iK 2 HZ . 


(C.103) 


Thus K 2 satisfies a closed equation and can be determined without knowing the 
unitary part U. This is particularly useful for quantifying the squeezing power of S, 
which only depends on K. 


C.9 Examples and Applications 


C.9.1 Single-Mode Case 


In the single-mode case, any Hamiltonian matrix can be written as 


H — uj q / T CcqX T co 2 Y , 


(C.104) 


where cu 0 , uq, uj 2 are real parameters, / is the identity matrix, and X and Y = iXZ 
are Pauli matrices. 11 From Eq. (C.37), we have the symplectic invariant 


tr (ZHZH) = (uj 0 I - oqX - cu 2 Y) (u 0 1 + uqX + w 2 Y) 


2 2 2 
UJ 2 


( 'l 'A 2 \ 2 

= sgn(o; 0 -uj 1 -uj 2 )uj , 


(C.105a) 

(C.105b) 

(C.105c) 


where u = \J\luq — \. Depending on the sign of this invariant, the canonical 

form of H bifurcates: the Hamiltonian H is equivalent to uicJa for ojq — — oj 2 > 0, 

or to iu(a 2 — a } 2 )/2 for Uq — uf — u 2 < 0. 

To take the Hamiltonian matrix H into canonical form, we introduce the symplectic 
matrix (Bogoliubov transformation), 


B = e 


iZM 


11 


Because Z* A XZX, the Pauli-Z is not included. 


(C.106) 
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Here the matrix M in the exponential shares the form of Eq. (C.104), and the whole 


generator thus takes the form 

iZM = 7 {k 3 X + k 2 Y — in 3 Z) , (C.107) 

where 7 > 0, and k 1 , k 2 , and k 3 are real numbers satisfying \k± + k 2 ~ K i\ — 1- Note 
that the following identity holds, 

(k]X + k 2 Y — in 3 Zy = (k\ + — k\) I = ±/ . (C.108) 

For + k 2 — k\ — 1, we have 
00 j 

B = (KiX + k 2 Y - m 3 Z) j (C.109a) 

5=0 J ' 

( °° 2 j \ / 00 2fc+l \ 

+ ( C ' 109b > 

= / coshy + (fiyX + k 2 Z — m 3 Z) sinliy . (C.109c) 

Similarly, for n\ + nl — k 3 = —1, we have 

B = I cosy + (k]X + k 2 Y — in 3 Z) siny . (C.110) 

We divide the symplectic transformation into two steps, 


B = B sqz (d)B rot (6), (C.lll) 

where B mt is a rotation and £> sqz a squeezer. The rotation part takes the form 

e„,(») = e'" z/2 , (C.112) 

which is both symplectic and unitary; it brings H into the form 

H‘ = BU») H B„ t (0) = w 0 / + + X , (C.113) 

where tan 0 = —oj 2 /uj 1 . The squeezer, which is both symplectic and Hermitian, takes 
the form 

B sqz (d) = . 


(C.114) 
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Applying B sqz to H', we have 

H'» = Bt qz (d) H' B sqz (d) = + V(#)X , (C.115) 

where £(0) = uj 0 and 7/(0) = yfu f + u j\. To determine the functions £($) and rj(d), 
we notice that H' a satisfies the differential equation 

d ^f = \ ( XH '» + H '»X) . (c.ne) 

which leads to the equations 

m = r](d) , m = m . (c.117) 

The solutions to Eqs. (C.117) depend on the initial conditions, and we discuss them 
separately. For coq — uj\ «*■ > 0, i.e., £(0) > 7/(0), we have 


£($) = a; cosh (0 + d ), rj(d) = a; sinh (0 + d ), 


(C.118) 


where tanh (p = y/wf+~tuf /oj q . For cug — uj 2 — 00% < 0, i.e., £(0) < 7 /( 0 ), we have 


£($) — uj sinh (0 + d ), rj(d) = uj cosh (0 + d ), 


(C.119) 


where tanh <p — u 0 / yj + uj\ . By choosing the squeeze parameter d = —0, we have 

the canonical form 


a 


uj I for ujq — uj\ — 0 J 2 > 0 , 
c uX for ujq — uj 2 — 0 J 2 < 0 . 


(C.120) 


The symplectic eigenvalues of H are real in the first case, and H is symplectic 
equivalent to the rotation Hamiltonian (phase shift); the symplectic eigenvalues are 
pure imaginary in the second case, and H is symplectic equivalent to the squeezing 
Hamiltonian. When ujq — uj 2 — CU2 = 0, it could be the trivial case H = 0, or 


H vn = I - X = 


'psq 




(C.121) 


where the subscript means p 2 for the reason we will see soon. Note that H psq is 
nilpotent, i.e., ( ZH psq ) 2 = 0, which means that it only has one eigenstate, (1, 1) T , 
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whose eigenvalue is zero; therefore, it is impossible to bring this Hamiltonian matrix 
into one of the canonical forms (C.85). The exponentiation of /7 psq , however, is quite 
easy because of the nilpotent condition, 

exp ( -itZH psq ) = / - itZH psq . (C.122) 

The evolution is linear in time instead of oscillatory or exponential, because the 
Hamiltonian describes a free particle, 

77 psq = - (a^a + acJ — a 2 — a^ 2 ) = p 2 . (C.123) 


C.9.2 Decomposing a Squeeze Operation into a Product of 
Rotations and the e~ ltp Operation 


/ 2 f 2 \ /q 

Sometimes, the squeeze operation <S(y) = e iya ~ a ' is not readily available for certain 
physical systems, but the rotation 7 Z{9) = e~ te ^ a a+aa and the e~ ltp operation are. 
In this subsection, I show how to decompose the squeeze operation into a product of 
the other two. 


The rotation, squeeze, and p 2 Hamiltonians take the forms 
77rot = - (a^a + aa , 

?4qz = ^{a 2 - a t2 ) , 

-U 1 / f | t 2 f 2 \ 2 

Tfpgq = -(a'a + aa 1 — a — a 1 ) = p , 


(C.124a) 

(C.124b) 

(C.124c) 


where we always use symmetric ordering of the annihilation and creation operators. 
The corresponding Hamiltonian matrices are 


1 0 


0 


Hm [o i) ' [i 0 ) ’ [-1 i)' (C,125) 

The evolution matrices, exponential in the Hamiltonian matrices, take the forms 

•SrotW = e~ l6ZHlot = I cos 6 — iZ sin 9 , (C.126a) 

■S sq z(7) = e~ nZHsqz = I coshy — X sinhy , (C.126b) 

S psq (t) = e ~ itZH = / - itZ /7 psq = / -itZ-tY . (C.126c) 
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Multiplying S rot and S psq together, we have 

SrotW 5 psq (t) 

= (/ cos 6 — iZ sin #)(/ — itZ — tY) (C.127a) 

= /(cos# — /sin#) — iZ(tcos9 + sin#) — tY cos 9 + tX sin# . (C.127b) 

By choosing t = — tan#, # G (— 7 t/ 2 , 7 t/ 2 ), we have 

5 rot (#) S psq (— tan #) = / sec # + (Y cos 9 — X sin #) tan # . (C.128) 

The terms with X and Y can be combined by inserting two rotations of opposite 


angles, 

5 r ot(-#/2 + 7t/4) S rot (#) S psq (-tan#) S rot (#/2 - 7T/4) 

= e i(0/2 -^ /4)z (/ sec # + (V cos # - X sin #) tan #) e -*#V 2 -V 4 ) z (C. 129a) 

= / sec# — Xtan# . (C.129b) 

By choosing tan# = sinh 7 , we have 

■5 sqz (7) = I coshy — X sinh 7 (C.130a) 

= 5 rot (#/2 + 7t/4) S psq (- tan#) S rot (#/2 - 7t/4) , (C.130b) 


Thus, we have decomposed the squeeze operator into three terms, each of which is 
either a rotation or an e~ ltv operation. 

C.9.3 Decoupling the Two-Mode Squeeze Hamiltonian 

Here, I show how to use the symplectic transformation to decouple the following 
two-mode squeeze Hamiltonian into a sum of two single-mode squeeze Hamiltonians, 


K tms = + b^b + 1) + iy(ab — a)b *). 


(C.131) 
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In matrix form, we have 



^ oj 0 0 — 


(a\ 

- (a} b t a b] 

0 oj —iy 0 


b 

a^ 

2 V / 

0 iy oj 0 



\iy 0 0 oj ) 


77 


and the corresponding Hamiltonian matrix is 


H 


tms 




(C.132) 


(C.133) 


where / and X are the 2x2 identity and Pauli-X matrices. Consider the following 
symplectic transformation of a beamsplitter, 


B bs P = exp 



(x 

o\ 


1 / 

-4^1 

A 

x ) 


~V2 ( 



(C.134) 


which brings the Hamiltonian matrix to 


t 1 I + iX 

dT u d _ _ 1 

D bsp n tms D bsp 2 y ^ 


0 


ojI -i 7 X \ I-iX 


I-iXj \d/X ul 


0 


0 

I+ iX 


(C.135a) 


( uil 7 / 

7 / uil 


(C.135b) 


The Hamiltonian matrix (C.135b) represents two identical independent single-mode 
squeezers. Although this problem can be done without any knowledge of symplectic 
methods, when the Hamiltonian matrix becomes more complicated, e.g., the frequen¬ 
cies oj for the two modes are different, it is quite difficult to simplify just by looking 
at the Hamiltonian, and the symplectic methods come to be a useful tool. 


C.9.4 Symplectic Eigenvalues for Two Modes 

In this subsection, I discuss how to solve the symplectic eigenvalue problem of the 
two-mode case (n = 2) without solving the symplectic eigenvalue equation (C.53). 
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It turns out that the two symplectic invariants, Eqs. (C.37) and (C.38), are 
particularly useful, 


Xd 

det (H ), 

Xt = 2 tr (ZHZH ). 

(C.136) 

Denoting the four eigenvalues of ZH by {±A 1; 

±A 2 }, we have 


|A!l 2 |A 2 | 2 

, |A x | 2 + A 2 2 , 

X 1 and A 2 are both real, 


-|A!| 2 |a 2 | 2 

, |Ail 2 -|A 2 | 2 , 

X 1 is real and A 2 pure imaginary, 

Xd; Xt < 

|AiI 2 |a 2 | 2 

, -(|Ar| 2 + |A 2 | 2 ), 

X 1 and A 2 are both pure imaginary, 


A 4 , 

2( Re A 2 - Im A 2 ), 

Ai = A 2 = A are complex. 


(C.137) 

Note that we leave out the case in which is pure imaginary and A 2 real, because 
can be obtained from the second case by swapping X 1 and A 2 . Assuming all the 
eigenvalues are nonzero, we can tell the type of the eigenvalues by the following 
procedure: (i) if y d < 0, H has a pair of real and a pair of pure imaginary eigenvalues; 
(ii) if Xd > 0 and Xt < 4y'd) H has a quadruplet of complex eigenvalues; (iii) if y d > 0 
and Xt > 4Xd and Xt > 0, H has two pairs of real eigenvalues; (iv) if Xd > 0 and 
Xt > 4y d an d Xt < 0, H has two pairs of pure imaginary eigenvalues. 

After determining the classes of the eigenvalues, we can solve for them by using 
Eq. (C.137). Knowing the symplectic eigenvalues, however, is not equivalent to 
knowing the canonical forms (C.85); for example, both ±H real have the same set of 
symplectic eigenvalues. 

C.9.5 Positiveness of the Hamiltonian Matrix 

In some cases, one wants to know if a quadratic Hamiltonian is stable (has a well 
defined ground state). For that purpose, one need only check whether the Hamiltonian 
matrix is positive definite [see Eq. (C.87)]. One way to do that is by the following 
condition for matrices of the spinor form (C.15): 

H is nonnegative if and only if > 0 and H 2 = \[fi\ C (y/Jfi) * 
for some contraction C. 12 


(C.138) 




Appendix C. Symplectic Methods for Quantum Mechanics 


147 


This theorem is pretty handy, because it reduces the positivity problem of a 2 n x 2 n 
matrix to that of an n x n matrix. 


C.9.6 Fluctuations Allowed by Quantum Mechanics 


For applications like quantum metrology or quantum computation, we may want to 
suppress quantum fluctuations to prevent noise. Zero quantum fluctuation is not 
possible because of the Heisenberg uncertainty principle (the covariance matrix cannot 
vanish in quantum mechanics). A natural question to ask is the following: what are 
the constrains on a given positive definite matrix R, so that there exist a quantum 
state whose covariance matrix is /?? This is an easy question for the single-mode case. 
The determinant of the given positive definite matrix must be greater than half of the 
Planck constant h: for the multi-mode case, however, the answer is not so obvious. 


In the basis of the quadratures, the covariance matrix of a quantum state is a 
2 n x 2 n real matrix, 


R = 


R\ Rl 

Ro Rn 


(C.139) 


where R 1Jk = 2 (AxjAxQ, R 2Jk = (Ax k Apj + ApjAxQ, and R 3Jk = 2 (ApjApQ 
with Axj = Xj — (Xj ) and Ap ] = p- — (p). In the basis of the annihilation and 


creation operators, we have 


A a,- = —-= (Ax,- 

3 V2 


iAp.j ), 

and the covariance matrix takes the new form 


A at = (Ax,- — iAp 


3 ) ’ 


(C.140) 


c = 




Ri R t 2 
R 



(C.141) 


where C 1 j k = (AajAa^, + Aa\Aaj), C 2 j k = 2 (AcifAa],). Note that C\ = C, and 
C 2 = C 2 , and the covariance matrix C has the same symmetry as the Hamiltonian 
matrix H. The good news is that the covariance matrix C is always positive definite, 


,t 


Ha A' 




12 


A contraction is an operator whose operator norm is less than one. 
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therefore, all of its eigenvalues are real. Thus, it can be brought into the following 
Williamson normal form by a symplectic transformation, 


= B ] CB = 


diag(A 1 ,..., \ n ) 0 

0 diag(A 1; ..., A n ) 


(C.142) 


Quantum mechanics constrains C can to satisfy 

A j = (AbjAbj + AbjAbj ) > 1, for j — 1, 2,..., n , 


(C.143) 


where Abj and Aare the annihilation and creation operators in the new basis, with 
their mean removed. On the other hand, we can always find a quantum state whose 
variance matrix is Eq. (C.142) if all its diagonal elements are greater than 1. Thus, 
given a 2 n x 2 n matrix C which satisfies O = C and C* = XCX , we have 


C is the covariance matrix for some quantum state if and only if 
its canonical form C can is greater or to than the identity /. 


(C.144) 


This condition, which requires the canonical form, is not convenient to use; an 
equivalent one is 

C can >Z. (C.145) 

Because Z is symplectically invariant, we have 

C > Z . (C.146) 

Note that this condition Eq. (C.146) does not require the canonical form, and it 
automatically excludes the pure-imaginary-eigenvalue and complex-eigenvalue cases. 

For pure Gaussian states, the equality in Eq. (C.143) is satisfied, i.e., 

C can = 6CB t = / <=*► C = ZB ] BZ (C.147) 

Consequently, we have 

C~ L = B ] B = ZCZ, (C.148) 


which implies the condition (C.146). 
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For the two-mode case, if we already know C > 0, then it is easy to restate the 
condition Eq. (C.144) in terms of symplectic invariants [see Subsec. C.9.4 for details], 

Xd — det(C) — A 1 A 2 , x,= ^tr(ZCZC) = A? + A|, (C.149) 

where A x > A 2 > 0 represent the positive eigenvalues. Since X 1 and A 2 must be greater 
or equal than one, we have 

At ~ \J Xt ~ 4y d = Af + A 2 — | Af — A 2 | > 2 . (C.150) 

This inequality, a necessary and sufficient condition for C to be physical, can be 
simplified to 

At<Xd + l- (C.151) 

For more detailed discussion on the two-mode case, see [PSL09]. 


C.9.7 Dynamics of the Covariance Matrix 


Consider a quantum state evolving under the quadratic Hamiltonian (C.16), 


n 


IQ s)( HM m) )(*) 

2 K > \H 2 (t) Hi (t) J y a y 


(C.152) 


The time evolution of the annihilation and creation operators in the Heisenberg 


picture is 



where the evolution matrix 5 satisfies 


(C.153) 


S(t) = - iZH(t)S(t ) . (C.154) 

The covariance matrix (C.141) at time t thus takes the form 

C(t) = S(t)C{0)S\t) , (C.155) 


and the governing equation is 

C(t) = -iZH{t)C(t)+iC(t)H(t)Z . 


(C.156) 
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Appendix D 

The Gross-Pitaevskii Ground 
State Is Stable 


A method is more important than a discovery, since the right method will 
lead to new and even more important discoveries. 

- Lev Landau 

In this Appendix, we investigate whether the ground state of the Gross-Pitaevskii 
Equation (GPE) is stable and show that it is because all the Bogoliubov excitations 
increase the energy of the condensate. The proof is constructed by noticing that both 
the perturbations of the condensate wavefunction and the Bogoliubov excitations are 
described by the Bogoliubov-de Gennes (BdG) equations [Bog58, dG66]. 

Consider a Bose-Einstein Condensate (BEC) govern by the Hamiltonian 

U = I ^ t (x)(-^-V 2 + H(x)jrl;(x) + | [rl; t (x)] 2 i|) 2 (x)hx . (D.l) 

The mean-field energy for the product state | ^(x) ) 0N is 

£„,= /iVV’*(x)(-Av 2 + V'(x))V’(x) + |]V(JV-l)[V(x)] 2 V’ 2 (x)dx. (D.2) 

The Gross-Pitaevskii (GP) ground state 0(x) is the state that minimize the mean-field 
energy E mi , and a perturbation around 0(x) takes the form 


0(x) -f y 1 - e 2 + e<p L (x) , 


(D.3) 
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where e is a small parameter, and is some normalized state orthogonal to 0(x). 

We neglect the change of an overall phase in Eq. (D.3), because it does not affect the 
mean-held energy. The first order perturbation of the mean-held energy for the GP 
ground state is zero, i.e., 5E^} = 0, which gives rise to the GPE 


(-Tv 2 + V(x) + g(N - l)«x)| 2 - k) 0(x) = 0 , 


(D.4) 


where /i is the chemical potential, and hereafter we neglect the difference between 
N — 1 and N. In the second order perturbation, the value of E mi cannot decrease for 
the GP ground state; i.e., the following quantity is nonnegative, 


SE, 


( 2 ) 

mf 


6 N 


II 


2 m 
12 


= / 01 (x) V“ + V(x) W> ± (x) - (j)*(x) + V(x) )0(x) 


r 


2 m 


2 + - 
2 


+ 9^ (2 |0(x) | |0 ± l 
With the GPE (D.4), we notice that 


(x)] 2 <#'i(x) +C.C.) - p(x) 


(D.5) 


dx 


^ = 1^1 0l(x)0i(x) dx — I f (x) (-Tv 2 + p( x ) + 9 ]V|^(x)| 2 )^(x) dx . 


(D.6) 


Putting Eq. (D.6) into Eq. (D.5), we have 


SE. 


( 2 ) 

mf 


= / 0i(x) + ^(x) - /^)^(x) + 2gN |0(x)f |0 ± (x)| 


e 2 N 


(D.7) 


+ ([0*( x )f 0±( x ) + c - c -) 

In matrix form, we have 

A 




( 2 ) 

mf 


e 2 N 


( 0_l I (01 


H^ + gNQWQ 


gNQtfQ* 


\ 


V 

where Q — 1 — | 0) (0 |, and 


gNQ*(<f>*) 2 Q H + gNQ*\(f){ z Q 


2 /"'i* 


/ 


|W| 

v l A )y 


(D.8) 


ft 2 


= + l " + 9 JV W -4*' 


(D.9) 


The GP ground state being stable is equivalent to the following matrix being positive 
definite in the subspace orthogonal to the condensate wavefunction | 0), 


^ neb 


2 * 


#gp + gNQmQ gNQ(f) z Q 


\ 


gN Q*(cp*) 2 Q H + gNQ*\4>\' z Q 


2 * 


(DIO) 


/ 
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Because the vector ( \(p±) |0l)) 7 in Eq. (D.8) is of a restricted form, an 

immediate conclusion that H nch is positive definite is not available. We notice, 
however, the following symmetry, 


H* nch = XH nch X , where X 




(D.ll) 


Suppose | q ) is an eigenvector of the Hermitian matrix H nchl we have 


H nch \q) = X\q) =► XH ncb XX\ q) — XX\q) (D.12a) 

=► H: ch X\q) = \X\q) (D.12b) 

=► H nch X\ q y = \X\ q y. (D.12c) 

Thus, both | q ) and X\ q)* are eigenvectors of H nch . More generally, any linear combi¬ 
nation of them is also an eigenvector of H ncb for the same eigenvalue A. Particularly, 
we introduce 


\q+) = \q) +X\q)* and \q_) = i(\q)-X\q)*) , (D.13) 

which satisfy the condition 


<Z±) =X\q ± Y ■ 


(D.14) 


Using this construction, we can always find a complete set of eigenvectors of H nch 
satisfying Eq. (D.14), and this is the same restriction on the vector ( \<j>±) | 4>*l ) ) T . 

Thus, Eq. (D.8) says that the minimum eigenvalue of H Ilch is nonnegative; i.e., H nch 
is positive definite. 


We can ask another question: Is the GP ground state the only GP eigenstate that 
is stable under Bogoliubov excitations? The answer is NO, as one simple example 
shows. Consider a BEC confined in a one-dimensional interval, x € [0, L], with 
periodic boundary condition. Any momentum eigenstate, i.e., ipk{ x ) = e kx /\[L. is 
also an eigenstate of the GP equation, 

/ h? < 9 2 2 \ 

\ 2mdx 2 + gN ^ )} ^ ~ ^ k )^ k = 0 ’ where = + gn ’ ( D - 15 ) 
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with n = N/L. Stability requires that 


H = 

r7 ncb 


Ho- p + gn Q k gn Q k e 2lkx Q ^ 


gnQ* k e 2lkx Q k H +gnQ* k 


(D.16) 


be nonnegative, where Q k (Q* k ) is the projector onto the orthogonal space of | ) 

(| V’-fc ))• Because the matrix H ncb only couples the vector ( | c+p ) | ^ k+p ) ) T /a/2 

to ( | t/’fc-p ) | V’fc-p ) ) T /v / 2 for p ± 0, we have the following in the subspace spanned 

by these two vectors, 

\ 

——-— + an an 

lip) _ 

h 2 (p 2 -2kp) 


U\P) — 

n neb ~ 


/ h 2 (p 2 + 2kp ) 

V 


gn 


2m 


+ 


(D. 17) 


'in I 


which is positive definite if and only if 


h 2 p 2 /2m + gn > yj ( h 2 kp/m ) 2 + gn h 2 k 2 /m < h 2 p 2 /Im + gn . (D.18) 


Stability requires that this condition holds for all p, and the strongest of which is 
when L —>■ oo, 


\k\ < yjgnmjh . (D.19) 

Since the right hand side of Eq. (D.19) is of finite value, a BEC can be stable even if 
it is not in the GP ground state. 
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Appendix E 

The Autonne-Takagi Factorization 


The mathematical facts worthy of being studied are those which, by their 
analogy with other facts, are capable of leading us to the knowledge of a 
physical law. 

- Henri Poincare 

In this Appendix, I discuss the Autonne-Takagi factorization theorem (see Corollary 
4.4.4 (c) in [HJ13]) in a way that physicists might found more approachable. The 
Autonne-Takagi factorization is useful for writing down the Schmidt decomposi¬ 
tion form of an arbitrary two-boson wavefunction [PY01]. Also, it is the tool to 
decomposing a multimode squeeze operator into single-mode squeezers (see Sec. C.8). 

The Autonne-Takagi factorization theorem says that any complex symmetric 
matrix $ = $ can be factorized into the following form: 

U t $U = A, (E.l) 

where U is a unitary matrix and A = diag(A 1; A 2 ,..., A n ). Note that the phases of 
the XjS are arbitrary; one choice is to make the AjS nonnegative by absorbing half of 
these phases into the column vectors of U and half to the row vectors of U . Consider 
the following time-reversal eigenvalue equation, 

$ | Uj) = A j | Uj )* , 


(E.2) 
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where the column vector | Uj ) is normalized, and j Uj)* = \ u*) is the complex 
conjugate of | Uj ) in the computational basis | ej ). The eigenvalue Xj = Xj | e 1 is 
generally complex, but the phase e l6j can always be absorbed into the eigenvector, 

$(e~ ldj/ ' 2 | Uj )) = |Aj-| (e ~ l6j/2 1 Uj ))* ; (E.3) 

without loss of generality, we assume A^ > 0 for j = 1,2,..., n. The transpose of the 
eigenvalue equation (E.2) reads 

(u* k |$ = X k (u k | , (E.4) 

and we have the following condition for the eigenvectors | Uj ) and | u k ), 

(u* k |$ | Uj) = Xj(u* k | u*) = X k (u k | Uj) . (E.5) 

For A j 7 ^ X k the above condition is equivalent to 

{u k \uj) = 0, (E. 6 ) 

and thus {| Uj) \ j = 1 , 2 ,... n} form an orthonormal basis when there are no 
degenerate eigenvalues . 1 With the eigenvectors, the matrix takes the form 

* = < E - 7 > 

j 

which is manifestly symmetric. Defining the unitary matrix U — Jjp | u 3 ) ( e. ; |, we 
have 


U T W = y I e 3 *) (^ I = £ VI e,) ( e } | = A , (E.8) 

3 3 

where the fact | e* ) = | ) is used. 

One interesting question is what equation determines the eigenvalues X y By using 
Eq. (E.2) twice, we have 

Uj ) = $*\j\u*j) = \Xj\ 2 \uj ) . (E.9) 

X The generalization to the degenerate case can be done by linearly combining the 
eigenvectors of the same eigenvalue. 
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Thus, | Uj ) is an eigenvector of the nonnegative matrix and we have the following 
equation that determines its eigenvalues, 

det(<f )t $ — |A| 2 1) = 0 , (E.10) 

where 1 is the identity matrix. 

To end this Appendix, I discuss how to validate the eigenvalue equation (E.2) 
from the conventional eigenvalue equation 

&<h\u) = = |A| 2 |-u). (E.ll) 

By introduce the vector | v) — $* | u *), Eq. (E.ll) takes the form 

$|u) = |u*), and <I>*| v* ) = |A| 2 | u) . (E.12) 

Note that | v ) is also an eigenvector of with the same eigenvalue | A| 2 , 

&$\v) = = |A| 2 r|u*) = |A| 2 |u) . (E.13) 

As a consequence, | v) is a multiple of | u) if the matrix is nondegenerate. Note 
also that we have the following normalization relation for the two vectors, 

<u|u) = (v*\v*) = (mI^Im) = \\\ 2 {u\u) . (E.14) 

Thus we can always choose the phase of A such that | v) = A| u ), and Eq. (E.12) is 
then equivalent to Eq. (E.2). 
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